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

Подгорнова Ольга Владимировна

ПОСТРОЕНИЕ ДИСКРЕТНЫХ ПРОЗРАЧНЫХ

ГРАНИЧНЫХ УСЛОВИЙ ДЛЯ АНИЗОТРОПНЫХ И

НЕОДНОРОДНЫХ СРЕД

Специальность 01.01.07 – вычислительная математика

АВТОРЕФЕРАТ

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

Москва – 2008

Работа выполнена в отделе №4 асимптотического и численного анализа задач механики сплошной среды Института прикладной математики им. М.В. Келдыша РАН.

Научный руководитель:

доктор физико-математических наук Софронов Иван Львович

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

доктор физико-математических наук, профессор, Рябенький Виктор Соломонович кандидат физико-математических наук, Холодов Ярослав Александрович

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

Объединенный институт геологии, геофизики и минералогии имени А.А. Трофимука СО РАН

Защита состоится «16» октября 2008 г. в 1140 часов на заседании диссертационного совета Д 212.156.05 при Московском физико-техническом институте (государственном университете), по адресу: 141700, Московская область, г. Долгопрудный, Институтский пер., д. 9, ауд. 903 КПМ.

С диссертацией можно ознакомиться в библиотеке МФТИ(ГУ).

Автореферат разослан «15» сентября 2008 г.

Ученый секретарь диссертационного совета кандидат физико-математических наук Федько О.С.

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

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

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

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

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

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

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

Поглощающие условия или ABC (исторически сложившееся название от английского Absorbing Boundary Conditions), относящиеся к локальным НИГУ, и использующие на границе дифференциальные операторы высокого порядка, требуют небольших вычислительных ресурсов, но не всегда устойчивы, и также не обладают достаточной точностью.

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

Прозрачные (или точные) граничные условия – ПГУ – относятся к классу нелокальных НИГУ и обеспечивают как высокую точность, так и устойчивость вычислений на большие времена. Они основаны на точных представлениях решений исходных уравнений в отбрасываемой внешней области и потому безупречны с математической точки зрения.

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

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

Еще одним результатом работы является аналитическое решение задачи построения ПГУ для широко используемой модели анизотропии – вертикально поперечно-изотропной среды (или VTI от английского Vertical Transverse Isotropy). В связи с этим отметим, что развитие аналитических ПГУ сдерживается относительной сложностью математического аппарата, из-за чего ПГУ остаются в нише фундаментальных исследований и распространены пока что далеко не на все возможные приложения.

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

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

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

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

• метод численно исследован на двух модельных задачах; продемонстрированы ожидаемые точность, устойчивость и экономность граничных условий при расчетах на длительные времена.

Развитый в работе метод впервые позволил распространить концепцию ПГУ на анизотропные и неоднородные среды.

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

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

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

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

На защиту выносятся следующие положения:

1. Метод построения приближенных дискретных прозрачных граничных условий для анизотропных и неоднородных сред.

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

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

Публикации Результаты исследований по теме диссертации изложены в восьми печатных работах, в том числе трех [2,3,4] из списка изданий, рекомендованных ВАК РФ.

В работах с соавторами лично соискателем выполнено следующее:

[1,2,5] – окончательно сформулирован метод, разработан и реализован алгоритм построения дискретных прозрачных граничных условий, проведены численные эксперименты на двух модельных задачах; [3,6] – предложен и реализован способ конструкции дискретных прозрачных граничных условий, непривязанный к численному методу решения основной задачи;

[7] – скорректирован алгоритм построения аппроксимации дискретных граничных условий, получены результаты лучшей точности; [8] – проведено сравнение спектрального и конечно-разностного подходов для нахождения дискретной функции Грина внешней задачи; [4] – реализован оператор прозрачных граничных условий для линеаризованной системы уравнений Эйлера, проведены численные эксперименты.

Апробация Результаты, полученные в работе, докладывались на конференциях:

• международной конференции «Workshop on nonlinear approximations in numerical analysis» (Москва, 2003), • XV и XVI Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики с приложением к многопроцессорным системам», посвященная памяти К.И. Бабенко (Абрау-Дюрсо, 2004, 2006), • На международном семинаре «Days on diffraction» (Санкт-Петербург, • на Всероссийской научно-практической конференции «Вычислительный эксперимент в аэроакустике» (Светлогорск, 2006), • на Всероссийской конференции по вычислительной математике «КВМ-2007» (Новосибирск, 2007), • на международной конференции «Matrix methods and operator equations» (Москва, 2007).

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

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

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

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

Во второй главе формулируется метод построения дискретных прозрачных граничных условий (ДПГУ), описываются численные алгоритмы их реализации, и проводится численное исследование условий на двух модельных задачах.

В параграфе 2.1 описываются модельные задачи: первая (A) – для анизотропной однородной среды, вторая (B) – для неоднородной изотропной среды. Обе задачи решаются в полярной системе координат ( r, ), а граница, для которой в дальнейшем строятся ДПГУ, является окружностью r = R.

A. Рассматривается процесс распространения акустических волн в движущейся среде. Анизотропия, проявляющаяся зависимостью скорости звука от направления распространения, вызвана наличием выделенного направления движения среды. Процесс описывается так называемым волновым уравнением в движущейся среде, получающимся из обычного волнового уравнения Здесь c – скорость распространения волн в покоящейся среде, 0 a < c скорость движения среды, обе величины ( c и a ) являются постоянными. В полярных координатах уравнение (2.1) имеет более сложный вид, в частности, содержит смешанные производные по времени и пространству.

Процесс распространения акустических волн в двухслойB.

ной среде описывается волновым уравнением с переменным коэффициентом (отвечающим за скорость распространения волн в В дальнейшем задачи A, B записываются единообразно, а именно где оператор c,a обозначает оператор Лапласа в движущейся или слоистой среде соответственно. Предполагается, что начальные данные u0, u1 и правая часть g имеют носитель, находящийся внутри области r R.

В параграфе 2.2 вводится понятие точных ДПГУ и описывается формальный алгоритм их построения.

Работа ведется с дискретным аналогом задачи (2.3), построенным для некоторого, определенного заранее способа дискретизации уравнений (разностная схема) на фиксированных сетках по времени и пространству ( и h соответственно). Обе сетки предполагаются равномерными с шагами количество точек сетки h по переменной. Для дискретизированных операторов tt и c,a вводятся обозначения Dtt и,,a, дискретное решение записывается как u,h, а его значения в точке t = t p, r = rj, = m обозначается через u jp,m. Таким образом, дискретный аналог дифференциальной задачи (2.3) имеет вид:

Точные ДПГУ определяются как условия, заданные на границе r = R и обеспечивающие совпадение решения дискретной задачи (2.4) в неограниченной области и этой же дискретной задачи, рассмотренной внутри Доказывается, что точные ДПГУ можно сконструировать из решений набора дискретных внешних начально-краевых задач, имеющих вид:

Здесь p – это дискретная дельта функция по времени, 0 = 1, p = 0 при p > 0 ; ( ), n = 0,1,…, M 1 – это некоторый полный дискретный базис на границе r = R (базис из синусов и косинусов обычно является удовлетворительным). Стоит отметить, что решения задач (2.5) есть не что иное, как проекции на базис дискретной функции Грина внешней начально-краевой задачи Дирихле. Показывается, что точные ДПГУ определяются соотношением между значениями решения при r = R и при r = R + hr. Действие ДПГУ заключается в разложении решения на границе r = R по базису ( ), дальнейшей свертки каждого элемента разложения с соответствуюn щей «базисной» дискретной функцией Грина и обратной сборке решения по базису:

Соотношение (2.6) используется для обновления решения на вспомогательной границе r = R + hr при наличии найденного по разностной схеме для внутренней области ( r R ) решения при r = R. В формуле (2.6) использованы обозначения F, F 1 для прямого и обратного разложения по базису ( ), символ обозначает дискретную свертку по времени, а матn рица G определяется через коэффициенты разложения G ( ) по базису ( ) разом где элементы G зависят от временного слоя t p.

В параграфе 2.3 описывается способ аппроксимации ДПГУ по пространству для уменьшения количества сверток в формуле (2.6).

Предлагается аппроксимировать матрицу G (которая в случае a является полностью заполненной) урезанной матрицей G с элементами G n,k при n nTBC и n k nband, при этом все остальные элементы полагаются равными нулю (см. Рис. 1). Параметры nTBC и nband подбираются на практике экспериментально.

гладкости решения, что позволяет считать вклад высоких G = гармоник по угловой переменM nTBC ной небольшим. Исключение n k nband предполагает малость этих элементов, наблю- Рис. 1. Схематичное изображение аппрокдающуюся при умеренной ани- симации матрицы G по пространству.

зотропии (для изотропии эти элементы равны нулю, поскольку матрица G диагональна).

В параграфе 2.4 описывается аппроксимация оператора ДПГУ по времени, основанная на алгоритме, предложенным A. Arnold, M. Ehrhardt, I.

Sofronov (2003) для построения дискретных граничных условий для уравнения Шредингера.

Каждый элемент матрицы G (2.7) аппроксимируется по времени суммой экспонент:

понент, Ln,k – количество экспонент (квадратные скобки используются для обозначения возведения в степень).

Использование аппроксимации (2.8) позволяет вычислять свертку в ДПГУ (2.6) по рекуррентным формулам В отличие от прямого вычисления свертки, где на каждом шаге по времени приходится перерабатывать всю историю решения на границе r = R, формулы (2.9) оперируют только со значениями на текущем шаге по времени и дополнительными «накопительными» переменными I ln,k. Такая численная локализация нелокальных ДПГУ существенно сокращает объемы вычислений и является ключевой составляющей метода.

Приближение (2.8) находится при помощи аппроксимации Паде, в которой используются только первые 2 Ln,k коэффициентов последовательноp сти G n,k. Следует отметить, что алгоритм Arnold, M. Ehrhardt, I.

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

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

(б) Абсолютная величина погрешности аппроксимации: (A) при использовании всех получаемых алгоритмом аппроксимации экспонент; (B) при исключении растущих экспонент.

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

Для анализа погрешности ДПГУ вычисляются три решения:

• решение U DTBC, полученное с использование ДПГУ на r = R ;

• решение U EXT, полученное на той же сетке что и U DTBC, но в расширенной области, так что отражения от внешней границы не достигают области r R за рассматриваемое время моделирования;

• эталонное решение U REF, вычисляемое в расширенной области на очень подробной сетке.

Решение U DTBC содержит как ошибку аппроксимации разностной схемы, так и ошибку, вызванную неточностью ДПГУ, а U EXT содержит только ошибку аппроксимации; решение U REF считается точным.

Ошибкой граничных условий называется величина, обусловленная только погрешностью аппроксимации матрицы G :

Совместная ошибка схемы и граничных условий есть Ошибка аппроксимации схемы определяется как Оператор ДПГУ не ухудшает точность аппроксимации разностной схемы, если максимальная (по времени) ошибка U appr + DTBC ( t ) близка к ошибке U appr. В частности, если U DTBC существенно меньше погрешности аппроксимации разностной схемы (например, более чем на порядок), то условие U appr + DTBC U appr выполнено.

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

В параграфе 2.6 оцениваются вычислительные затраты, необходимые для реализации ДПГУ.

Ресурсы, необходимые для вычисления и аппроксимации матрицы G, велики, однако эти вычисления достаточно провести только один раз (сетка, разностная схема, внешняя граница r = R, параметры среды для внешней области фиксированы), а затем использовать для расчетов множества внутренних задач. При этом допускается изменение любых параметров внутренней задачи: источников, начальных данных и даже уравнений и их параметров, при условии, что изменение сосредоточено строго внутри При использовании точных ДПГУ количество операций растет линейно с увеличением номера временного слоя p и квадратично зависит от количества точек M по переменной, что, конечно, неудовлетворительно.

В то же время, в ДПГУ, аппроксимированных по предлагаемому в диссертации алгоритму, количество операций не изменяется с номером p, линейно зависит от количества экспонент L = max Ln,k, зависимость от M также обычно линейна.

В параграфе 2.7 приводятся результаты численного исследования ДПГУ для волнового уравнения в движущейся среде.

Рассматривается задача с нулевыми начальными данными и возмущением, порожденным гладкой правой частью. Проводится набор экспериментов с использованием центрально-разностной схемы второго порядка точности при фиксированной скорости c = 1 и положении внешней границы R = 2 для двух значений параметра a, a = 0.2 и a = 0.7 на трех вложенных сетках G1, G 2, G3, каждая из которых в два раза мельче предыдущей по обоим направлениям.

Численно проверено, что L2 -нормы элементов G n,k убывают при удалении от диагонали G (см. Рис. 3), что подтверждает правомерность используемой пространственной аппроксимации ДПГУ (см. параграф 2.3).

Для демонстрации точности ДПГУ проводятся расчеты для набора параметров nTBC и nband. Анализ результатов подтверждает, что при достаточно больших значениях nTBC и nband ДПГУ не ухудшают точность разностной схемы. При этом в рассматриваемых экспериментах для самой подробной сетки ( G3 ) точность разностной схемы 103, а для самой грубой сетки ( G1) точность порядка пары процентов. Точность ДПГУ U DTBC ограничена снизу погрешностью аппроксимации G по времени. Приводятся численные эксперименты, демонстрирующие, что точность аппроксимации улучшается при увеличении количества экспонент L.

Рис. 3. Нормы элементов матрицы G для сетки G1. (а) a = 0.2, (б) a = 0.7.

Предлагается способ апостериорной оценки точности временной аппроксимации, основанный на вычислении погрешности в частотной области. На примерах демонстрируется, что при вычислении ошибки необходимо исключать из рассмотрения высокие частоты (см. Рис. 4). Мотивацией является наблюдение, что высокие частоты, в которых сосредоточена основная ошибка при аппроксимации G n,k суммами экспонент, не присутствуют в точном решении.

В параграфе 2.8 приводятся численные эксперименты с использованием ДПГУ для волнового уравнения в слоистой среде. Анализ ДПГУ проводится аналогично уже рассмотренному случаю движущейся среды, демонстрируются результаты такого же порядка точности. За исключением некоторых деталей, в частности более медленного убывания норм G n,k при удалении от диагонали, поведение и свойства ДПГУ сохраняются.

В параграфе 2.9 обсуждаются результаты второй главы.

Рис. 4. Точность аппроксимации G ( L = 40, круглые значки) и G100 ( L = 100, треугольные значки), для каждого элемента G n,k как функция номера диагонали n k для сетки G 3. (а) L2 -норма по времени (совпадает с нормой во всем частотном диапазоне). (б) Точность после исключения высоких частот.

В третье главе предлагается и численно исследуется способ построения ДПГУ со сглаженной функцией Грина.

основная идея метода, поясняются его преимущества и недостатки. 0. главе ДПГУ зависят от разностной схемы, что обусловлено наличием 0. времени в краевых условиях задач симости предлагается заменить Рис. 5. Локальный сплайн t. Точками дискретную дельта функцию на ее гладкий аналог s ( t ). Требуется, чтобы с помощью s ( t ) можно было бы аппроксимировать гладкие функции с порядком не меньшим чем порядок разностной схемы относительно. Гладкость s ( t ) позволяет перейти от дискретной постановки к дифференциальной и тем самым обеспечить независимость граничных условий от разностной схемы (зависимость от параметров сетки hr и сохраняется). Дополнительно G, в силу гладкости, аппроксимируется меньшим коs личеством экспонент, что уменьшает расходы на реализацию граничных условий, однако приводит к возрастанию предварительных расходов на вычисление G s.

В параграфе 3.2 описываются детали метода построения ДПГУ со сглаженной функцией Грина. Построение проводится на сетках и h, определенных в главе 2. В качестве s ( t ) используется локальный сплайн В.С. Рябенького, построенный по значениям дискретной (см. Рис. 5).

Такая функция непрерывна вплоть до 3-ей производной, имеет локальный носитель [ 2,2 ] и обеспечивает аппроксимацию произвольной достаточно гладкой функции с порядком 3.

Сглаженная функция Грина G s строится аналогично дискретной G (2.7), при этом используются решения набора внешних начально-краевых задач:

Доказывается, что для гладкой функции f, являющейся суммой не более M гармоник ( задачи аппроксимируется с порядком 3 на слое r = r1 дискретной функцией u,h = F 1 G s F f,h. Делается вывод, что использование ДПГУ вида разностных схем.

Результаты численных экспериментов с использованием ДПГУ со сглаженной функцией Грина для волнового уравнения в движущейся среде приводятся в параграфе 3.3. Рассматривается только случай a = 0.7 и сетки G3.

Для нахождения сглаженной функции Грина задачи (3.1) решаются численно по разностной схеме, используемой в дальнейшем и для решения внутренней задачи (центральные разности второго порядка). Для корректного описания гладкой s ( t ) сетка G3 измельчается во всех направлениях в d раз. Элементы сглаженной G приведены на Рис. 6 для двух значений d, d = 4 и d = 8. Вычисленные Gsn,k являются гладкими функциями с незначительными осцилляциями, вызванными недостаточной точностью решения задач (3.1), при измельчении сетки (увеличении d ) амплитуда осцилляций уменьшается.

Продемонстрирована высокая точность ДПГУ со сглаженной функцией Грина, и уже d = 4 является достаточным для обеспечения точности ~ 103.

В параграфе 3.4 исследуется возможность сокращения предварительных расходов, требуемых на вычисление сглаженной функции Грина G s за счет использования высокоточной разностной схемы для решения внешних задач (3.1). В качестве альтернативы центральным разностям второго поряка точности рассматривается спектральный метод, использующий псевдоспектральную дискретизацию по пространству и аппроксимацию по времени, основанную на разложении в ряд по обобщенным функциям Лагерра (Конюх, Михайленко, Михайлов, 2001).

DGF DGF

10, Рис. 6. Абсолютная величина элементов Gs для n = 10 и k = 10 (а), k = 11 (б). Показаны элементы дискретной G (DGF) и гладкой матрицы G (MDGF), вычисленные на сетке, измельченной по отношению к исходной в d = 4 и d = 8 раз. Временной интервал соответствует 300 временным шагам.

Исследование проводится численно на модельном примере волнового уравнения, где внешняя двумерная задача распадается на набор одномерных задач для каждой Фурье гармоники m Сравнение проводится на примере одного параметра m ( m = 6 ).

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

SPECTRAL

ет, что для рассматриваемой зада- number of operations чи, несмотря на спектральную скорость сходимости, метод, использующий функции Лагерра и псев- пространству, является более ресурсоемким, чем обычная явная мация второго порядка вплоть до Рис. 7. Количество операций относительточностей 106 (см. Рис. 7). но точности решения. Сплошная кривая Четвертая глава посвящена соответствует конечно-разностному поданалитическим прозрачным гранич- ходу, пунктирная – спектральному.

ным условиям (ПГУ).

В параграфе 4.1 прозрачные граничные условия для линеаризованной системы уравнений Эйлера (Софронов И.Л., 1999) применяются в аэроакустической задаче моделирования звукопоглощающей конструкции.

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

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

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

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

Уравнения линейной упругости формулируются в виде системы первого порядка для вектора скорости и тензора напряжений, которая для осесимметричного случая описывается уравнениями движения:

и законом Гука, имеющим вид:

Здесь ( ur, u z ) и ( rr,, zz, rz ) – это компоненты вектора скорости и тензора напряжений соответственно, – плотность, A11, A12, A13, A33, A44 – параметры среды.

Рассматривается процесс распространения волн в бесконечной полуполосе, 0 r +, Z min z Z max, с граничными условиями:

на верхней и нижней границах ( z = Z min, z = Z max ). Вычислительной областью является ограниченный прямоугольник [ 0, R ] [ Z min, Z max ], на правой границе которого r = R строятся прозрачные граничные условия, обеспечивающие прохождение волн без отражений.

При построении ПГУ для системы (4.1) (4.2) используется подход, изложенный в работах Софронова (1992-2006), Alpert, Greengard, Hagstrom (2000), где сначала граничные условия формулируются в спектральном пространстве (преобразование Фурье вдоль границы и преобразование Лапласа по времени), а затем преобразовываются обратно в физическое пространство.

Результирующая формула граничных условий в спектральном пространстве имеет вид:

где ur, rz – это ur, rz компоненты после преобразования Фурье по переменной z и преобразования Лапласа по переменной t, а s и l – соответствующие двойственные переменные. Матрица где скалярные функции ± определены как Здесь K1 – функция Макдональда первого порядка, а параметры a0, a2, b1, c1, c1,,,,, зависят только от параметров среды:

Для построения вычислительно эффективных ПГУ матрица P представляется в виде:

матрицы Q0, Q1, Q2 зависят только от параметров среды, матрицы Q1, K ( s ) зависят дополнительно как от номера Фурье гармоники, так и от положения границы R. Элементы матрицы K ( s ) стремятся к 0 при s +.

После обратных преобразований Фурье и Лапласа условий (4.3) с использованием представления (4.4) ПГУ в физических переменных принимают вид:

Здесь F обозначает преобразование Фурье, функция K ( t ) является обратным преобразованием Лапласа от K ( s ), символ обозначает свертку по времени.

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

В проведенных численных экспериментах граничные условия (4.5) используются совместно со схемой J. Virieux (1986) на сдвинутых сетках для системы (4.1), (4.2). На границе вычислительной области r = R компоненты ur и rz обновляются с использованием дискретизированных уравнений (4.5), в то время как для обновления остальных компонент используются дискретизированные уравнения основной системы (4.1), (4.2).

Работоспособность и точность полученных граничных условий проверяется на наборе одномерных тестов для первых пятнадцати Фурье гармоник ( l = 0,…,14 ). Волновой процесс, содержащий обе квази P- и квази Sволны примерно одинаковой амплитуды, возбуждается источником, присутствующим во всех компонентах тензора напряжений. По времени источник описывается импульсом Риккера, по координате z определяется соответствующей Фурье гармоникой и является точечным по r координате. Рассматриваются две сетки: грубая – с 256 точками по r направлению и подробная – с 512 точками. Для оценки точности ПГУ используется Снорма на временном интервале моделирования, вычисленная около границы r = R.

Показывается, что точность расчетов с локальными граничными условиями невысока (достигает десятков процентов с увеличением номера гармоники, см. Рис. 9а), при этом она не улучшается при измельчении сетки.

Напротив, точность расчетов с использованием нелокальных ПГУ оказывается высокой (погрешность менее 2 103 для подробной сетки), при этом сеточная сходимость со вторым порядком схемы Virieux не нарушается (см. Рис. 9б).

Рис. 9. Относительная точность ПГУ для набора Фурье гармоник. (а) Локальные граничные условия. (б) Нелокальные граничные условия.

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

В приложении 1 приведены необходимые определения и свойства аппроксимации Паде.

В приложении 2 сведены необходимые данные об обобщенных функциях Лагерра.

В приложении 3 описан используемый в диссертации алгоритм разложения в ряд по функциям Лагерра.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

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

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

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

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

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

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

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Софронов И.Л., Подгорнова О. В. Нелокальные спектральные граничные условия для волнового уравнения в движущейся среде: Препринт / ИПМ им. М.В. Келдыша. – М. 2004. – №53. – 19 с.

2. Sofronov I.L., Podgornova O.V. A spectral approach for generating nonlocal boundary conditions for external wave problems in anisotropic media // J. of Scientific Comput. – 2006 – V.27, N 1–3. – P. 419-430.

3. Подгорнова О.В. Построение оператора дискретных неотражающих граничных условий для моделирования волн в движущейся среде // Математическое моделирование – 2007. – Т.19, № 8. – С. 75-82.

4. Александров А.В., Дородницын Л.В., Подгорнова О.В. Сравнение неотражающих граничных условий на примере задачи с внешним источником колебаний // Математическое моделирование – 2007. – Т.19, № 8. – С. 55-65.

5. Софронов И.Л., Подгорнова О.В. Спектральный метод построения слабоотражающих граничных условий для моделирования волн в анизотропных средах // Материалы XV Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики с приложением к многопроцессорным системам», посвященная памяти К.И. Бабенко – АбрауДюрсо, 2004. – С. 22-23.

6. Podgornova O.V., Sofronov I. L. Toward efficient numerical generation of low-reflecting boundary conditions for anisotropic media // Материалы международного семинара «Days on diffraction» – Санкт-Петербург – 2005. – P. 73-74.

7. Подгорнова О.В., Софронов И.Л. Эффективное вычисление и аппроксимация функций Грина для оператора слабо-отражающих граничных условий в анизотропной среде // Материалы XVI Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики с приложением к многопроцессорным системам», посвященная памяти К.И. Бабенко – Абрау-Дюрсо, 2006. – С. 19-20.

8. Podgornova O.V., Sofronov I.L. Laguerre-Pseudospectral Discretization vs. Finite-Differences for Discrete Green Function Calculation // Abstracts of the II International Conference «On matrix methods аnd operator equations» – Moscow, 2007. – P. 66-67.

ПОСТРОЕНИЕ ДИСКРЕТНЫХ ПРОЗРАЧНЫХ ГРАНИЧНЫХ

УСЛОВИЙ ДЛЯ АНИЗОТРОПНЫХ И НЕОДНОРОДНЫХ СРЕД