Компьютерное моделирование процесса динамического деформирования акустических, упругих и упруговязкопластических сред с упрочнением тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Ми Синь

• Ми Синь
• кандидат наук
• 2026

Оглавление диссертации кандидат наук Ми Синь

Введение

Глава 1. Механико-математические модели деформируемых сред

1.1 Акустическая среда

1.2 Линейно упругая изотропная среда

1.3 Случай криволинейной системы координат

1.4 Упруговязкопластическая среда с упрочнением

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

2.1 Сеточно-характеристический метод

2.1.1 Алгоритм для внутренних узлов

2.1.2 Реализация граничных условий

2.1.3 Реализация контактных условий

2.2 Метод многостадийного расщепления

2.2.1 Алгоритм расчёта во внутренних узлах

2.2.2 Алгоритм расчёта в граничных узлах

2.3 Криволинейные расчётные сетки

2.3.1 Расчётные сетки специального вида

2.4 Явно-неявная схема

2.4.1 Частичная неявная аппроксимация

2.4.2 Случай линейных зависимостей

2.5 Описание расчётной программы

Глава 3. Результаты моделирования

3.1 Эмпирическое исследование порядка сходимости

3.2 Многослойные геологические среды

3.3 Упрочняющиеся упруговязкопластические среды

Заключение

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

Стр.

Список рисунков

Список таблиц

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

Введение

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

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

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

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

В настоящее время на практике широко применяются как подходы, основанные на лучевом приближении, так и на решении волнового уравнения. Первая группа технологий включает такие подходы, как метод трассировки лучей и интегральный метод Кирхгофа, которые, в силу своей относительно низкой вычислительной сложности, получили широкое практическое распространение. Они подходят для приближённого решения задач миграции и инверсии. Теория лучей обеспечивает эффективное асимптотическое описание распространения волнового поля в высокочастотном пределе, её математические основы и физическая интерпретация подробно изложены в монографии Cerveny [1]. Вторая группа технологий основывается на полноволновых методах решения определяющей гиперболической системы уравнений. Основные уравнения распространения упругих волн и соответствующие им определяющие соотношения исходят из классической теории сплошной среды, а соответствующие положения можно проследить в работах Ландау и Лифшица по теории упругости [2] и Лойцянского по механике жидкости и газа [3]. В работе Igel [4] систематически рассматриваются методы численного моделирования сейсмических волн с вычислительной точки зрения, тогда как Chapman [5] проводит единый и глубокий анализ физических механизмов распространения сейсмических волн и соответствующих математических моделей. К применяемым на практике методам относятся: метод конечных разностей, метод конечных элементов, метод спектральных элементов. Они, несмотря на высокую вычислительную сложность, позволяют гораздо точнее моделировать процесс распространения сейсмических волн в неоднородных, анизотропных средах, при наличии сложных границ раздела слоёв.

Численное моделирование волновых процессов методом конечных разностей широко применяется в задачах сейсмики и механики сплошных сред. Исследования показали, что точность конечно-разностных схем и уровень численной дисперсии в значительной степени зависят от пространственного шага сетки и порядка аппроксимации схемы. Было установлено, что переход от схем второго порядка к схемам более высокого порядка аппроксимации позволяет существенно снизить требования к числу узлов, приходящихся на длину волны, и повысить точность моделирования волновых полей [6; 7]. Эти результаты послужили основой для дальнейшего развития конечно-разностных

методов высокого порядка аппроксимации и их анализа в рамках формальной теории устойчивости и сходимости, сформулированной в классических работах по численным методам математической физики [8]. С целью дальнейшего подавления численной дисперсии и повышения устойчивости были предложены различные оптимизированные схемы, включая методы с совместной дискретизацией по времени и пространству, а также неявные конечно-разностные схемы [9—11]. Оптимизация дисперсионных соотношений с использованием анализа распространения плоских волн и метода наименьших квадратов позволяет эффективно уменьшать дисперсию и использовать более крупные шаги по времени при сохранении приемлемой точности [9; 10]. Активно развиваются подходы, ориентированные на моделирование в сложных геометриях и неоднородных средах, включая многоблочные криволинейные сетки [12], методы переходного слоя для сопряжения жидких и твердых сред [13], а также методы погруженной границы на декартовых сетках [14], что существенно расширяет область применимости конечно-разностных схем в практических задачах сейсмического моделирования. При решении нелинейных гиперболических уравнений и систем важную роль играют высокоразрешающие схемы, обеспечивающие баланс между высокой точностью и подавлением не физических осцилляций вблизи разрывов решения. В этой связи были разработаны и исследованы схемы типа Годунова, центральные схемы, а также различные модификации WENO подходов [15—18]. Несмотря на их высокую эффективность при описании ударных волн и резких градиентов решения, вычислительная сложность таких методов часто ограничивает их прямое применение в задачах крупномасштабного волнового моделирования. Это обстоятельство стимулирует дальнейшие исследования, направленные на разработку высокоточных, устойчивых и вычислительно эффективных конечно-разностных схем, применимых к моделированию волновых процессов в сложных многокомпонентных моделях.

В исследованиях, связанных с применением метода конечных элементов (МКЭ) для задач волновой динамики, различные исследователи расширяли и совершенствовали различные аспекты: стратегии дискретизации, вычислительную эффективность, методы постановки искусственных граничных условий. В работе [19] были использованы нелинейные определяющие соотношения и проведено численное моделирование осесимметричных нагруженных вращающихся оболочек в рамках смешанной МКЭ задачи. В работе [20] представлен раз-

рывный МКЭ повышенного порядка для задачи трёхмерного распространения упругих волн. В ней используется пространственная дискретизация уравнений первого порядка в скоростях и напряжениях разрывным методом Галёркина и параллельные алгоритмы, что позволило эффективно моделировать поле упругих волн в сложных трёхмерных моделях. В работе [21] систематически проанализировали влияние пространственной и временной дискретизации на точность расчётов с использованием МКЭ, особенно для распространения высокочастотных волн в упругих и упруго-пластических средах, и предложили рекомендации по улучшению традиционных критериев дискретизации. Для повышения вычислительной эффективности моделирования больших областей и сейсмических волн в работе [22] предложен усовершенствованный метод МКЭ со значительным снижением потребления памяти и PML-граничным условием. В области развития реализации искусственных поглощающих границ Matzen [23] предложил свёрточную реализацию PML в МКЭ для динамических уравнений упругости в смещениях. В работе [24] была разработана реализация PML в МКЭ для эффективного моделирования волн в анизотропных средах.

Метод спектральных элементов (Spectral-Element Method, SEM) сочетает в себе возможность описания сложных геометрий и высокую точность спектрального элемента. В работе [25] описан вариант SEM на неструктурированных треугольных сетках с использованием полиномиальных аппроксимаций высокого порядка для точного моделирования распространения упругих волн в сложных геометрических областях (наличие рельефа), систематически обсуждаются проблемы диагонализации матрицы масс и вычислительной эффективности. В обзорной работе [26] представлены теоретические основы и численные реализации SEM для региональных и глобальных сейсмических симуляций, показана возможность генерации высокоточных синтетических сейсмограмм в трёхмерных неоднородных моделях Земли и выделен потенциал применения метода для инверсии и томографии. В работе [27] подробно описана полная реализация SEM для трёхмерного моделирования распространения упругих волн, с акцентом на структуру диагональной матрицы масс, полученной с использованием интегрирования по Гауссу-Лобатто-Лежандру, и на эффективные стратегии параллельных вычислений, что подтвердило стабильность и точность SEM при масштабных трёхмерных численных симуляциях.

При численном моделировании сейсмических волн выбор конкретной модели среды непосредственно влияет на получаемые характеристики колебаний и

достоверность проводимого моделирования. Акустическая модель среды - наиболее упрощенная модель, учитывает только плотность и модуль объёмного сжатия среды, не описывает сдвиговую деформацию и наличие сдвиговых волн. Она широко используется при проведении крупномасштабного предварительного моделирования и быстрого прямого расчёта, например, при моделировании шельфовой сейсморазведки. В моделировании и обработке сейсмических записей акустические модели широко применяются благодаря своей простоте и высокой вычислительной эффективности, однако их физическая интерпретация и область применения требуют критического анализа. В работе [28] аналитически исследованы особенности распространения волн в анизотропных акустических средах. Выявлено, что даже при нулевой скорости сдвига акустическое приближение может генерировать псевдосдвиговые волны с ненулевой групповой скоростью, осложняя использование методов сейсмической томографии. Эффективность акустических моделей в практическом сейсмическом картировании была продемонстрирована, например, в работе [29]. На основе данных DAS (Distributed Acoustic Sensing) для VSP (Vertical Seismic Profiling) реализована акустическая обратная временная миграция для сложных условий мерзлоты, обеспечивая высокое пространственное разрешение газогидратных резервуаров. Для повышения вычислительной эффективности моделирования в сильно неоднородных средах в работах [30—32] применили многомасштабные базисные функции МКЭ, что позволило снизить размерность задачи при сохранении информации о подсеточных масштабах. Для обратной сейсмической задачи был предложен метод Q-компенсации для акустической импедансной инверсии [33], позволяющий получать высокоточные оценки параметров среды.

В отличие от неё, линейно упругая модель дополнительно учитывает наличие у среды ненулевого модуля сдвига и может одновременно описывать как продольные волны (P-волны, волны сжатия-растяжения), так и поперечные волны (S-волны, волны сдвига). В современных программных комплексах по численному моделированию сейсмических волн она широко используется, особенно, при расчётах на высокопроизводительных вычислительных системах. В монографии [34] подробно изложены основные формы уравнений волнового движения в линейно-упругой среде, механизмы распространения объемных и поверхностных волн, математическое описание источника возмущений, а также характеристики волнового поля в слоистых и неоднородных средах применительно к сейсмологии. В обзорной работе [35] систематически рассматриваются

различия между идеальными упругими волнами и волнами в реальных материалах, с акцентом на граничные эффекты и нелинейные явления. Также перечислены основные приложения теории упругих волн в микросейсмологии, сейсмической разведке и инженерной сейсмике. С учётом неоднородности и нелинейности материалов в работе [36] выполнено численное моделирование распространения волн в слоистых нелинейно-упругих средах, используя высокоразрешающую конечно-объёмную схему, получено совпадение с экспериментальными наблюдениями. В работе [37] предложен метод моделирования поля упругих деформаций на основе метода точечного источника (базисного решения), систематически проанализированы ошибки и сходимость метода для различных краевых задач, продемонстрирован его потенциал для получения высокоточных решений линейно упругих задач. С точки зрения вычислительной реализации, в работе [38] построен численный алгоритм решения двумерного уравнения упругости на основе метода конечных разностей во временной области.

В последние годы высокоточное численное моделирование сейсмических волн [39—41] становится одним из перспективных направлений развития современных технологий сейсморазведки. Несмотря на значительные успехи, остаются не решённые вопросы в области минимизации ошибок пространственной дискретизации, снижения численной дисперсии, более точной аппроксимации граничных и контактных условий.

Для удовлетворения требований к высокой точности моделирования было проведено множество усовершенствований расчётных алгоритмов. Например, метод конечных разностей высокого порядка [42—44], разрывный метод Галер-кина [45—49] и метод многошагового расщепления операторов [50—52] нашли широкое применение, поскольку они позволяют эффективно подавлять численную дисперсию и ложные отражения (численные артефакты, немонотонности) при сохранении высокой вычислительной эффективности. В работе [42] для уравнений упругости второго порядка на криволинейной расчётной сетке построены высокоточные конечно-разностные операторы с узким шаблоном на основе метода ЗиттаМоп-Ву-Раг^Б (БВР). С помощью энергетических оценок систематически доказана устойчивость полученных разностных схем. В работе [43] предложена схема конечных разностей четвертого порядка как по пространству, так и по времени для уравнения упругости второго порядка с переменными коэффициентами. Численные эксперименты показали

значительное увеличение допустимого шага по времени и преимущество перед традиционными схемами второго порядка. В последние годы разрывный метод Галеркина (Discontinuous Galerkin, DG) благодаря сочетанию высокой точности и гибкости работы с геометрией в сложных и сильно неоднородных средах постепенно стал важным инструментом высокоточного численного моделирования сейсмических волн. В работе [45] была использована идея временной аппроксимации ADER (Arbitrary high-order DERivatives) в рамки DG, формируя метод ADER-DG для двумерных неструктурированных треугольных сеток, что позволило обеспечить согласованную высокую точность по пространству и времени для численного решения линейно упругих уравнений. Явное разрешение разрывов на границах элементов и введение численного потока на основе решения задачи Римана обеспечили сохранение высокой точности и устойчивости в сложных неоднородных средах. Далее метод был систематически расширен на трёхмерные неструктурированные тетраэдральные сетки для моделирования распространения упругих волн в изотропных средах, при этом численные эксперименты [46] показали сохранение спектральной сходимости и значительное снижение вычислительных затрат при заданной точности. На этой основе было проведено расширение метода ADER-DG для вязкоупругих сред с использованием модели обобщённого тела Максвелла, что позволило моделировать реальные эффекты затухания [47]. Для анизотропных сред в работе [48] было проведено обобщение метода ADER-DG, а численные эксперименты подтвердили высокую точность и устойчивость на неструктурированных сетках. В дальнейшем Dumbser [49] предложил сочетание p-адаптации с локальным шагом по времени (Local Time Stepping, LTS) в ADER-DG, что позволило использовать различные степени полиномов и временные шаги для разных элементов, значительно повысив вычислительную эффективность в условиях сильных изменений масштаба сетки и сложного трёхмерного рельефа, обеспечивая практическую применимость для крупных сейсмических моделей.

Современные технологии построения расчётных сеток, такие как криволинейные структурные сетки [53—59] и адаптивные сетки [60—62] значительно улучшили способность моделировать сложные геометрии геологических включений и разломов. В то же время, в противовес классическим численным методам, активно развиваются методы машинного обучения, включая физически информированные нейронные сети [63; 64], которые предоставили быстрые

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

В последние годы сеточно-характеристический метод (Grid Characteristic Method, GCM) [65—68], использующий математические особенности гиперболической системы уравнений, получил широкое развитие в области моделирования сейсмических волн. Его основы достаточно подробно были изложены в работе [69]. Его высокая вычислительная эффективность и адаптируемость к широкому кругу задач были продемонстрированы как при использовании структурированных [70—72], так и неструктурированных расчётных сеткок [73]. Для повышения порядка аппроксимации многомерных сеточно-характеристических схем возможно применение метода многошагового расщепления [52].

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

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

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

1. Построение сеточно-характеристического метода пятого порядка аппроксимации на криволинейных расчётных сетках для двумерных динамических задач акустики и упругости.

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

3. Развитие явно-неявных схем для динамических задач упруговязкопластических сред с эффектом упрочнения.

4. Разработка исследовательского программного комплекса для решения задач динамического деформирования рассмотренных сред.

Научная новизна:

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

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

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

4. Получены численные решения трёхмерных задач динамики упруговяз-копластических сред с упрочнением.

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

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

применяемых языках C++ и Python и протестированы на известных аналитических решениях.

Основные положения, выносимые на защиту:

1. Построены двумерные сеточно-характеристические схемы пятого порядка аппроксимации для расчёта процесса распространения акустических и упругих волн на прямоугольных и специального вида структурированных расчётных сетках.

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

3. Разработан программный комплекс для моделирования волновых эффектов в деформируемых средах. Получены численные решения двумерных и трёхмерных задач.

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

Апробация работы. Основные результаты работы докладывались на следующих конференциях:

1. «66-я Всероссийская научная конференция МФТИ», г. Долгопрудный,

2024 г.

2. «Марчуковские чтения - 2024», г. Новосибирск, 2024 г.

3. «Quasilinear Equations, Inverse Problems and their Applications - 2024», Сириус, 2024 г.

4. «67-я Всероссийская научная конференция МФТИ», г. Долгопрудный,

2025 г.

5. Научная конференция «Численное моделирование в механике сплошных сред», посвященная 100-летию академика О. М. Белоцерковского, г. Долгопрудный, 2025 г.

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

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

Публикации. Основные результаты по теме диссертации изложены в 4 печатных изданиях [75—78], в том числе 3 из них индексируются в Web of Science/Scopus: [75—77].

Объем и структура работы. Диссертация состоит из введения, 3 глав, заключения и 0 приложен. Полный объём диссертации составляет 93 страницы, включая 15 рисунков и 18 таблиц. Список литературы содержит 98 наименований.

Глава 1. Механико-математические модели деформируемых сред

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

Список литературы диссертационного исследования кандидат наук Ми Синь, 2026 год

// - 1))

: Б^1

+

+БП

^ (/Б^Гб™(кп)

/Бп : Бп

1

= 68Т+1,

(2.36)

кп+1 = кп + 1

Р2

((' (

/ (^+4

1

+

(' (

/ (кп)

1

. (2.37)

Сворачивая уравнения (2.36), (2.37) последовательно с б

ёП+1 сп тт ^-+1

ц и б:

получаем систему уравнений относительно сверток X = //БП+1 : Бп ,г = ф

. -^п+1 у =

: б

'^П+1 . оП+1

: Бе"- с вычисленными значениями

Т = : Бп =

х/Б^Б^1, 2 = у/БГ+1 : БГ+1.

Эта система приводит к формулам:

6Х + (X// (Г+1) - 1)) = 65,

(2.38)

в2^п+1 - (X// (Г+1) - 1)) = в2^п,

(2.39)

где ёп+1 = - ^{Е(г//Р-1)), кп = кп + ^ {Г (Т// (кп) - 1)), ¡3 = л/£п+г : ё^1, Т = /8п : 8п.

Легко видеть, что в этом случае промежуточный девиатор ёп+1 и его свертка ¡^ вычисляются на основе результатов "упругого"временного шага.

Из уравнений (2.36)-(2.39) получаем формулы для локальной корректировки значений напряжений в виде

6§П+1 цП+1

^ = 5 + )) = ~3rX, (2.40)

кп+1 = 6 -Х^) /в2 + кп. (2.41)

Финальные значения 8п+1 и кп+1 вычисляются как решение (2.38)-(2.39), подставленное в условия (2.40)-(2.41). Предельный переход 6 ^ 0 приводит к

8п+1 = ,п+1 { Г (Т/,1 (кп) - 1)) 8 =8е - Т 6 .

Для этого уравнение (2.38), справедливое на п + 1-м временном слое, переписывается для п-го временного слоя:

6Т + { Г (Т// (кп) - 1)) = 6/8пТ~8п.

Отсюда

{ Г (Т// (кп) - 1)) = 6 ^/8пТ§п - ,

а в формуле для промежуточного девиатора, учитывая 6/в2 = О"о/ (2ц), можно устранить особенность при 6 ^ 0:

8п+1 = 8«+1 - , ^ - : ,

кп+1 = кп + ^ -Х^~кп = кп + (/8п : 8п - /8п : .

Этот прием, предложенный в свое время в работах [86; 87], широко использовался в вычислительной практике, например, в [88—91].

2.4.2 Случай линейных зависимостей

Рассмотрим специальный случай системы (2.38)-(2.39), когда функция упрочнения имеет вид /(к) = 1 + ак, где параметр а ^ 0 регулирует режим и скорость изменения предела текучести. При а > 0 имеем увеличение предела текучести (упрочнение), при а = 0 имеем идеальную вязкопластичность. Пусть дополнительно Г(х) = х.

При этих предположениях система уравнений (2.38)-(2.41) имеет вид:

6Х + /-—Х~~г - Л = 65, \ 1 + акп+1 /

Р2кп+1 - (Т+^ТТ - ^ = в>к"

8п+1 / „ \

8п+1 = ^Х, кп+1 = 6 (¡3 -X) /в2 + кп. ¡к

Исключая из этой системы кп+1, получаем уравнение для X: X = (^1 + 6 (§ -Х^ + ва6 (¡к - ^ + ак.

Для удобства решения введем новую переменную С = £> - X и новую нотацию для комбинации параметров у = а/в2. Уравнение для С примет вид:

у62С2 + {1 + 6 {1+ У + аГ)) С -(¡3 - 1 - аГ) = 0.

Аналитическое решение с С > 0 равно:

^{1 + 6 (1+ У + акп))2 + 4у62 - 1 - акп^ - (1 + 6 (1+ у + ак;п^ ^ = 2^6 .

Следовательно,

^{1 + 6 (1+ У + аЬ))2 + 4у62 - 1 - аЬ) - (1 + 6 (1+ у + акп)) X = ^ 2^6

Формулы для девиатора и параметра упрочнения на следующем временном слое примут вид:

/

sn+1 =

gn+l

S-

+ ak^j

2y62

V

(1 + 6(1 + Y +

2Y62

Ln+1 = A в2

To

^ ^1 + 6(1+ Y + + 4y62 - akn^j

2y62

V

(1 + 6(1+ Y +

2Y62

Эти формулы для случая малого коэффициента упрочнения а и малого

параметра Y ^ 1 преобразуются к виду:

1 + акп + 6 ¡3(1+ Y +

X

sn+1

1 + 6 (1 + y +

gn+l 1 + akkn + 6S? (1 + Y +

S ' '

(2.40)

1 + 6 (1 + y +

kn+1

ao .S - 1 - akn 1 + 6 (1 + y + ok"

+ kn.

Они допускают предельный переход к чистому упругопластическому материалу:

g

Zn+1

S

(1 + аГ) ,

kn+1 « (S - 1 - + kn.

Условие реализации режима упрочняющейся УВП: S - 1 - akn ^ 0.

2.5 Описание расчётной программы

Исследовательский программный комплекс Яее^ развиваемый на кафедре информатики и вычислительной математики МФТИ, представляет собой крупномасштабную платформу для численного решения гиперболических систем дифференциальных уравнений в частных производных, предназначенную для моделирования процесса распространения волн в неоднородных деформируемых средах. Программное обеспечение имеет модульную архитектуру, интегрируя в себе множество физических моделей, таких как акустика, линейная упругость, трещиноватые и флюидонасыщенные среды. Оно оснащено продвинутым механизмом генерации расчётной сетки и высокопроизводительным вычислительным ядром, опирающимся на параллельные алгоритмы ОрепМР/МР1/СиЭЛ, обеспечивая единую и эффективную вычислительную среду для сложных вычислительных экспериментов.

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

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

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

3. Для акустической и изотропной линейно упругой моделей реализованы двумерные сеточно-характеристические схемы пятого порядка аппроксимации на специального вида криволинейных структурных сетках. Для них расчёт Якобиана координатного преобразования проводится

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

4. Реализован программный алгоритм многошагового операторного расщепления 5-го порядка аппроксимации по времени. Для каждого пространственного направления X,Y, задаётся s = 7 коэффициентов шага по времени а. После последовательного применения одномерных сеточно-характеристических операторов получаем численное решение на промежуточном временном слое. После повторения процедуры s раз итоговое численное решение соответствует временному слою п + 1.

5. Реализован программный модуль, осуществляющий упруговязкопла-стическую корректировку с учётом эффекта упрочнения. На первом шаге алгоритма с использованием сеточно-характеристического метода строится решение изотропной линейно упругой задачи. Затем проводится шаг корректировки значений тензора напряжений в тех узлах, в которых выполнен критерий выхода из упругого режима:

S - 1 - акп ^ 0.

Все указанные функциональные расширения были выполнены на языке программирования C++ самостоятельно автором диссертации, включая этапы тестирования и проведения последующих вычислительных экспериментов с использованием доработанного программного комплекса Rect и скриптов на Python.

Глава 3. Результаты моделирования

3.1 Эмпирическое исследование порядка сходимости

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

Рассматривалась однородная акустическая среда размерами 1000 х 1000 м. Плотность р задавалась равной 1000 кг/м3, скорость акустической волны ср — 1500 м/с. На боковых и нижней границах расчётной области использовалось неотражающее граничное условие. На всей верхней границе области прикладывалось заданное давление в виде

р(х,у = 0,г) = вт6(2п/г),г ^ 0,

с / = 10 Гц. Решение этой задачи в момент времени Т = 0,5 е можно найти аналитически как

V 750 < ж < 1000, Р (х,у,0.5)={ п

в г п (75х - 20п), 0 < 750

Были проведены расчёты с последовательным измельчением расчетной сетки. Различие между аналитическим и численным решениями в значениях Р измерялось в нормах Ь1 и

ИЛк =^2 IЛI • ||/||ьто =тах I ¡31.

3

Численный порядок сходимости а определяется таким образом:

, ||Е2^1 Ьр .

а = 1 од2 ,, ,, ^, р=1, ж. ^РьуЬр

где Ь - шаг по пространству, Е2н - ошибка численного решения при размере шага сетки 2Н, Е^ - ошибка численного решения при размере шага сетки Н.

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

Таблица 2 — Анализ сходимости численного решения задачи с граничным условием в однородной акустической среде.

Mesh size, h Error in L\ Error in Ьж Order in L\ Order in Lж

norm norm norm norm

10.0000

5.0000

2.5000

1.2500

0.6250

1.1660E+05

4.8804E+04 4.5292E+03 1.7514E+02 5.5672E+00

4.4912E-01 1.7514E-01 2.5643E-02 1.0064E-03 3.1990E-05

1.257 3.430 4.693 4.975

1.359 2.772 4.671 4.975

Был исследован порядок сходимости построенной двумерной схемы на задачах о распространении наклонных под углом а плоских Р- и Б-волн в акустической и линейно-упругой средах (см. Рисунок 3.1). Размер расчётной области составлял 600 х 600 м. Физические характеристики для акустики задаются р = 1000 кг/м3 и скоростью акустической волны ср = 2000 м/с, для упругости р = 1000 кг/м3, ср = 2000 м/с и с3 = 1000 м/с. Время расчета Т = 0.5 с. Результаты, представленные в Таблицах 3 и 5 получены с использованием стандартного расщепления и схемы пятого порядка. В свою очередь, результаты, приведённые в Таблицах 4 и 6 получены на основе разработанного в данной работе многошагового расщепления в сочетании со схемой пятого порядка. Зависимости ошибки численного решения от мелкости разбиения расчётной сетки при фиксированном числе Куранта приведены на Рисунке 3.2 (акустика, упругость).

Далее численно решалась двухслойная задача для подтверждения достижения повышенного порядка сходимости при наличии нормально падающей на контактную границу плоской волны. Размер всей расчётной области составлял 100 х 200 м. Первый слой задавался с плотностю ра = 1000 кг/м3 и скоростью акустической волны са = 1600 м/с. Второй: рь = 1200 кг/м3 и съ = 2000 м/с. Первоначально задавалась плоская волна на расстоянии 60 м от контактной границы, направленная вдоль оси ОУ. В этом случае прошедшая и отраженная волны могут быть рассчитаны аналитически. Коэффициент прохождения

Таблица 3 — Анализ сходимости одношаговой численной схемы для задачи с начальной Р-волной под углом а = 30° в акустической среде.

Mesh size, h Error in L1 Error in Ьж Order in Li Order in

norm norm norm norm

2.0000 1.0000 0.5000 0.2500 0.1250

4.1479E+08 7.5796E+07 1.3038E+07 3.2198E+06 8.0700E+05

4.9130E+05 1.1231E+05 1.9125E+04 4.5794E+03 1.1492E+03

2.452 2.539 2.018 1.996

2.219 2.554 2.062 1.995

Таблица 4 — Анализ сходимости многошаговой численной схемы для задачи с начальной Р-волной под углом а = 30° в акустической среде.

Mesh size, h Error in Li Error in Order in Li Order in

norm norm norm norm

2.0000 1.0000 0.5000 0.2500 0.1250

5.1314E+08 1.3443E+08 8.6838E+06 2.9673E+05 9.4008E+03

6.0974E+05 1.6667E+05 1.0911E+04 3.6512E+02 1.1515E+01

1.932 3.952 4.871 4.980

1.871 3.933 4.901 4.987

Таблица 5 — Анализ сходимости одношаговой численной схемы для задачи с начальной Я-волной под углом а = 30° в упругой среде.

Mesh size, h Error in Li Error in Order in Li Order in

norm norm norm norm

2.0000 3.6947E+08 4.2607E+05 - -

1.0000 6.1750E+07 8.5265E+04 2.581 2.321

0.5000 8.1832E+06 1.3532E+04 2.916 2.656

0.2500 3.2900E+06 4.8032E+03 1.315 1.494

0.1250 1.5380E+06 2.1386E+03 1.097 1.167

Таблица 6 — Анализ сходимости многошаговой численной схемы для задачи с начальной Б-волной под углом а = 30° в упругой среде.

Mesh size, h Error in Li Error in Lж Order in Li Order in Lж

norm norm norm norm

2.0000 1.0000 0.5000 0.2500 0.1250

4.4901E+08 1.1781E+08 7.6193E+06 2.6044E+05 8.2513E+03

5.3406E+05 1.4607E+05 9.5724E+03 3.2043E+02 1.0106E+01

1.930 3.951 4.871 4.980

1.870 3.932 4.901 4.987

Ят и коэффициент отражения Яд равны, соответственно,

2гь

Rt =

rr =

Za + Zb Zb — Za

^а +

Результаты проведённых тестов сходимости подтверждают сохранение схемой пятого порядка сходимости (см. Таблицы 7 и 8).

Таблица 7 — Анализ сходимости численного решения задачи при реализации контактного условия на границе раздела акустических сред с различными физическими параметрами для прошедшей волны.

Mesh size, h

Error in Li norm

Error in L norm

оо

Order in Li norm

Order in L norm

оо

1.000 7.5703E+08 6.2916E+05 - -

0.500 2.3430E+08 1.9814E+05 1.692 1.667

0.250 1.6211E+07 1.5632E+04 3.853 3.664

0.125 5.5549E+05 5.2805E+02 4.867 4.888

0.0625 1.7489E+04 1.6622E+01 4.989 4.989

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

Таблица 8 — Анализ сходимости численного решения задачи при реализации контактного условия на границе раздела акустических сред с различными физическими параметрами для отраженной волны.

Mesh size, h Error in Li Error in Order in Li Order in

norm norm norm norm

1.000 0.500 0.250 0.125 0.0625

1.1631E+08 3.9616E+07 3.0844E+06 1.0862E+05 3.4266E+03

1.1888E+05 4.0870E+04 3.5638E+03 1.2311E+02 3.8798E+00

1.554 3.683 4.828 4.986

1.540 3.520 4.855 4.988

Расчётная область покрывала прямоугольник размером 100 х 200 м. Деформированное состояние двухслойной среды описывалось системой уравнений изотропной линейной упругости с коэффициентами: ра = 1000 кг/м3, = 1600 м/с, ¿а) = 1100 м/с, ръ = 1200 кг/м3, с(рЬ) = 2000 м/с и с^ = 1200 м/с.

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

2 2(Ь) 11т = г

z\a) + z\b)

п = 2 - 2

= , - тип волны.

Последовательное измельчение расчётной сетки со множителем 1 позволило эмпирически оценить порядки сходимости схемы по нормам Ь1 и (см. Таблицы 9 - 12). Полученные данные подтвердили сохранение схемой пятого порядка сходимости для данных случаев.

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

Таблица 9 — Анализ сходимости численного решения задачи с начальной Р-волной при реализации контактного условия на границе раздела упругих сред с различными физическими параметрами для прошедшей волны.

Mesh size, h Error in Li Error in Lж Order in Li Order in Lж

norm norm norm norm

1.000 0.500 0.250 0.125 0.0625

7.5703E+08 2.3430E+08 1.6211E+07 5.5549E+05 1.7489E+04

6.2916E+05 1.9814E+05 1.5632E+04 5.2805E+02 1.6622E+01

1.692 3.853 4.867 4.989

1.667 3.664 4.888 4.989

Таблица 10 — Анализ сходимости численного решения задачи с начальной Р-волной при реализации контактного условия на границе раздела упругих сред с различными физическими параметрами для отраженной волны.

Mesh size, h Error in Li Error in Lж Order in Li Order in Lж

norm norm norm norm

1.000 1.1631E+08 1.1888E+05 - -

0.500 3.9616E+07 4.0870E+04 1.554 1.540

0.250 3.0844E+06 3.5638E+03 3.683 3.520

0.125 1.0862E+05 1.2311E+02 4.828 4.855

0.0625 3.4266E+03 3.8798E+00 4.986 4.988

Таблица 11 — Анализ сходимости численного решения задачи с начальной Б-волной при реализации контактного условия на границе раздела упругих сред с различными физическими параметрами для прошедшей волны.

Mesh size, h Error in L\ Error in Lж Order in Li Order in Lж

norm norm norm norm

1.000 0.500 0.250 0.125 0.0625

4.4723E+08 1.3952E+08 9.9399E+06 3.4195E+05 1.0787E+04

4.2514E+05 1.3550E+05 1.0589E+04 3.6487E+02 1.1488E+01

1.681 3.811 4.861 4.986

1.650 3.678 4.859 4.989

Таблица 12 — Анализ сходимости численного решения задачи с начальной Б-волной при реализации контактного условия на границе раздела упругих сред с различными физическими параметрами для отраженной волны.

Mesh size, h

Error in L

i

Error in L

оо

Order in Li

Order in L

оо

norm norm norm norm

1.000 4.9422E+07 4.9488E+04 - -

0.500 1.6050E+07 1.5769E+04 1.623 1.650

0.250 1.1494E+06 1.3258E+03 3.804 3.572

0.125 3.9857E+04 4.5164E+01 4.850 4.876

0.0625 1.2566E+03 1.4231E+00 4.987 4.988

представляет изучение более общего случая, когда направление распространения волнового фронта не совпадает с криволинейными линиями сетки. Было проведено исследование скорости сходимости численного решения задачи о распространении акустической волны в среде с плотностью р = 1000 кг/м3 и скоростью распространения продольной волны ср = 2000 м/с. Расчётная область покрывалась сеткой, получаемой дискретизацией данного аналитического отображения

х = £,,

у = п + ао£2.

(3.1)

где параметр а0 = 0.0005. Отметим, что заданием а0 = 0 можно получить исходную квадратную расчётную сетку размера 600 х 600 м.

Кроме того, были проведены вычислительные эксперименты с более общей моделью деформируемого твёрдого тела - изотропной линейно упругой моделью. В этом случае дополнительно использовалось значение параметра с3 = 1000 м/с. В обеих постановках фронт исходной (акустической/упругой) Р-волны отстоял на расстояние 400 м от нижней границы расчётной области и распространялся вертикально вниз (не вдоль линий сетки). Нормы ошибки рассчитывались через фиксированное число шагов по времени (см. Таблицы 13 и 14).

Таблица 13 — Анализ сходимости численной схемы для задачи с начальной Р-волной на криволинейной сетке в акустической среде.

Mesh size, h

Error in Li norm

Error in L norm

оо

Order in Li norm

Order in L norm

оо

2.0000 2.1456E+09 7.7196E+05 - -

1.0000 7.9590E+08 2.7508E+05 1.431 1.489

0.5000 6.7155E+07 2.6181E+04 3.567 3.393

0.2500 2.4542E+06 9.3466E+02 4.774 4.808

0.1250 7.8102E+04 2.9622E+01 4.974 4.980

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

Таблица 14 — Анализ сходимости численной схемы для задачи с начальной Р-волной на криволинейной сетке в упругой среде.

Mesh size, h Error in Li Error in Order in Li Order in

norm norm norm norm

2.0000 1.0743E+09 3.8619E+05 - -

1.0000 3.9865E+08 1.3765E+05 1.430 1.488

0.5000 3.3660E+07 1.3109E+04 3.566 3.392

0.2500 1.2305E+06 4.6806E+02 4.774 4.808

0.1250 3.9162E+04 1.4833E+01 4.974 4.980

расчётной области при использовании специального класса криволинейных расчётных сеток и семистадийной схемы расщепления (см. Таблицы 15 и 16).

Таблица 15 — Анализ сходимости численной схемы для задачи с начальной Р-волной, заданной под углом а = —5 на криволинейной сетке в акустической среде.

Mesh size, h Error in Li Error in Order in Li Order in

norm norm norm norm

2.0000 1.0000 0.5000 0.2500 0.1250

1.7507E+09 6.4658E+08 5.3794E+07 1.9591E+06 6.2373E+04

7.6346E+05 2.6924E+05 2.5158E+04 8.9504E+02 2.8357E+01

1.437 3.587 4.779 4.973

1.504 3.420 4.813 4.980

3.2 Многослойные геологические среды

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

Таблица 16 — Анализ сходимости численной схемы для задачи с начальной Р-волной, заданной под углом а = —5 на криволинейной сетке в упругой среде.

Mesh size, h Error in Li Error in Order in Li Order in

norm norm norm norm

2.0000 1.0000 0.5000 0.2500 0.1250

8.8399E+08 3.2676E+08 2.7199E+07 9.9093E+05 3.1550E+04

3.8501E+05 1.3581E+05 1.2706E+04 4.5211E+02 1.4325E+01

1.436 3.587 4.779 4.973

1.503 3.418 4.813 4.980

На первом этапе рассматривалась задача распространения акустической волны от точечного источника в многослойной акустической среде. Акустическая модель была построена на основе информации об упругой модели из статьи [92] путём удаления из параметров величны cs. Моделируемая область содержит восемь слоев с различной плотностью, скоростью распространения акустической волны и толщиной (данные приведены в Таблице 17). Горизонтальный размер модели составлял 4000 м. Источник был помещен в центр дневной поверхности. На границе задавалось давление в виде нескольких периодов периодического импульса

Р(t) = sin6 (2П).

Модель была покрыта квадратной расчётной сеткой с шагом 10 м. Временной шаг 40 мкс был выбран для удовлетворения условия устойчивости Куранта. На нижней и боковых границах области использовались неотражающие граничные условия. Пространственные распределения давления, полученные в момент времени Т = 780 мс двумя различными численными схемами, представлены на Рисунке 3.3. Он демонстрирует более высокое пространственное разрешение построенной численной схемы относительно часто используемой на практике схемы с интерполянтом третьей степени без многошагового расщепления [93; 94]. Все типы акустических волн - исходная сферическая, прошедшие, отражённые на контактных границах - чётко фиксируются в численном решении.

Таблица 17 — Акустические характеристики слоёв горизонтально-слоистой геологической модели.

Слой No. Плотность, Скорость Толщина, м

кг/м3 P-волны,

м/с

1 2000 2170 500

2 2300 2130 100

3 2200 2500 300

4 2300 2680 100

5 2400 3000 400

6 2700 5550 100

7 2800 6000 150

8 2850 6000 2000

На втором этапе рассматривалась задача распространения сейсмических волн от точечного источника в многослойной изотропной линейно упругой (каждый слой) среде.

Данные о модулях упругости и геометрических размерах каждого из 8 слоёв соответствовали информации из работы [92]. Они приведены в явном виде в Таблице 18. В центре верхней части расчётной области прикладывалось периодическое граничное условие, соответствующее внешней вертикально направленной силе со временной зависимостью вида:

f(x,y,t) = - sin6 (2nt),

до t ^ 50 мс.

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

Волновая картина напряжений иуу после Т = 780 мс физического времени представлена на Рисунке 3.4. Подтверждено чёткое получение объёмных продольных и поперченых волн, включая переотражённые сейсмические отклики от горизонтальных границ раздела слоёв. Отметим также отсутствие не физических осцилляций, что делает построенную схему перспективной для проведения прикладных компьютерных расчётов.

Таблица 18 — Упругие характеристики слоёв горизонтально-слоистой геологической модели.

Слой N0. Плотность, Скорость Скорость Толщина, м

кг/м3 Р-волны, Я-волны,

м/с м/с

1 2000 2170 674 500

2 2300 2130 795 100

3 2200 2500 1090 300

4 2300 2680 1220 100

5 2400 3000 1385 400

6 2700 5550 3144 100

7 2800 6000 1250 150

8 2850 6000 1550 2000

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

Было проведено численное моделирование процесса распространения сейсмических волн в трёхслойной геологической модели. Вся расчётная область размерами 90 х 150 м исходно состояла из трёх прямоугольных расчётных сеток. Далее была произведена их деформация специальным образом. Так, средняя подобласть полностью трансформировалась координатным преобразованием (3.1) с фиксированным значением а0 = 0.0004. Верхняя и нижняя подобласти также были трансформированы подобным преобразованием, однако, в них значения а0 постепенно уменьшались от а0 = 0.0004 до а0 = 0 при удалении от контактной границы. Таким образом, самая верхняя и самая нижняя границы расчётной области оставались полностью горизонтальными (см. Рисунок 3.5).

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

р = 1000 кг/м3, ср = 2000 м/с. В упругом случае дополнительно задавалось с, = 1000 м/с.

Была численно решена задача о вертикальном распространении продольной волны сжатия-растяжения, исходно заданной на глубине 20 м. Для выполнения условия устойчивости использовались расчётные параметры: = 100 мкс, Ах = 0.5 м. Рассчитывался процесс на протяжении 500 шагов по времени, что достаточно для прохождения волной обеих криволинейных контактных границ. Отметим, что тогда как в центральной расчётной сетке расчёт ведётся с пятым порядком аппроксимации, в двух других сетках порядок схемы снижается из-за зависимости параметра а0 от вертикальной координаты.

Сравнение волновых картин, полученных по различным расчётным схемам [95], представлено на Рисунках 3.6 и 3.7. Их анализ подтверждает более полное сохранение амплитуды распространяющейся волны, а также отсутствие значимых численных артефактов.

3.3 Упрочняющиеся упруговязкопластические среды

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

На первом этапе рассматривался случай а = 0, при котором среда описывается в рамках вязкопластической модели. Рассматривалась расчётная область размерами 50 м х50 м х10 км. Механические параметры задавались значениями р = 2500 кг/м3, ср = 4500 м/с, с3 = 2250 м/с, а = 112500 Па. Для обеспечения достаточной точности расчёта и устойчивости схемы для упругой части [85; 96—98] шаг пространственной сетки задавался равным 5 м, шаг временной сетки - 1 мс.

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

°зз = —---.

В частности, из

ап = ^22 = Л/(Л + 2^)азз,

получаем:

(ап - а/3)2 + (а22 - а/3)2 + (азз - а/3)2 = а2, а = ап + а22 + азз = 3г+г^азз.

Л + 2ц,

Окончательно,

Л + 2^

4

Для иллюстрации возможности расчёта процесса деформирования с заданным параметром т = 0 в работе была проведена серия тестовых расчётов. В них в широком диапазоне изменялось значение 6. При этом рассматривался только случай Г (х) = х. Полученные распределения полей напряжений представлены на Рисунке 3.8. Прямое сопоставление с аналитическим решением для среды без вязкости подтверждает корректность численного расчёта.

Далее исследовалось влияние величины скорости упрочнения упруговяз-копластической среды на процесс динамического деформирования. Параметр а варьировался в широких пределах, однако условие у << 1 всегда выполнялось. Результаты проведённых численных расчётов представлены на Рисунке 3.9. Видно, что при больших значениях 6 изменение параметра а влияет на скорость распространения волны пластической деформации.

Рисунок 3.1 — Начальное распределение волнового поля (наклонная волна): Р в акустической модели (сверху); ахх в упругой модели (снизу). Контуром показана область оценки численного решения.

Jb ^ТЪ —(L5 ChO 05

log h

To ^lb 00 (X5

log h

Рисунок 3.2 — Зависимость ошибки численного решения от мелкости разбиения расчётной сетки при фиксированном числе Куранта в log-log формате для акустической (сверху) и линейно упругой (снизу) среды.

Рисунок 3.3 — Волновое поле Р в многослойной модели: одномерная схема третьего порядка и одношаговое пространственное расщепление (сверху) и построенная в работе схема (снизу).

Рисунок 3.4 — Волновое поле иуу в многослойной модели: одномерная схема третьего порядка и одношаговое пространственное расщепление (сверху) и построенная в работе схема (снизу).

Рисунок 3.5 — Трёхслойная среда, покрытая криволинейными сетками.

10 20 30 40 50 60 70 80 90

Рисунок 3.6 — Поле давления Р в трёхслойной модели: одномерная схема третьего порядка и одношаговое пространственное расщепление (сверху) и построенная в работе схема (снизу).

О 10 20 30 40 50 60 70 80 90

100- 100 — 0.0е+00

- -20000

90 90 - -40000 - -60000

80- -80 - -80000

70- -70 - -100000 - -120000

60- -60 50 --140000 - -160000 —-180000

50- - -200000

40- -40 - -220000 - -240000

30- -30 - -260000 - -280000 х

20- -20 1- -300000 к -320000

10- -10 и -340000

|- -360000

0- - -380000

-10- -10 - -400000 - -420000

-20- -20 --440000 --460000

-30- -30 - -480000 - -500000

-40- -40 - -520000 - -540000

-50 -50 - -5.6е + 05

0 10 20 30 40 50 60 70 80 90

0 10 20 30 40 50 60 70 80 90

Рисунок 3.7 — Поле ихх в трёхслойной модели: одномерная схема третьего порядка и одношаговое пространственное расщепление (сверху) и построенная в

работе схема (снизу).

о.о

-0.5 -1.0

1Л

1/1 т г-

aj -1.5

с

о

2 -2.0

CD

£

b -2.5 -3.0

---Linear viscosity with i ---Linear viscosity with < ..... Linear viscosity with < - Analytical solution

5 = 10 5 = 10"2