Построение дискретных прозрачных граничных условий для анизотропных и неоднородных сред тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

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

□03447372

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

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

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

я О СЕН 2008

Москва-2008

003447372

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

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

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

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

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

кандидат физико-математических наук, Холодов Ярослав Александрович

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

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

Защита состоится ок-Г-^Г-8 2008 г. в ^ ^ часов

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

С диссертацией можно ознакомиться в библиотеке МФТИ(ГУ). Автореферат разослан« » С-е^-я. 2008 г.

Ученый секретарь диссертационного совета /?

кандидат физико-математических наук Федько О.С.

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

Актуальность темы

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

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

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

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

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

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

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

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

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

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

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

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

Еще одним результатом работы является аналитическое решение задачи построения ПГУ для широко используемой модели анизотропии -вертикально поперечно-изотропной среды (или 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» (Санкт-Петербург, 2005),

• на Всероссийской научно-практической конференции «Вычислительный эксперимент в аэроакустике» (Светлогорск, 2006),

• на Всероссийской конференции по вычислительной математике «КВМ-2007» (Новосибирск, 2007),

• на международной конференции «Matrix methods and operator equations» (Москва, 2007).

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

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

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

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

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

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

А. Рассматривается процесс распространения акустических волн в движущейся среде. Анизотропия, проявляющаяся зависимостью скорости звука от направления распространения, вызвана наличием выделенного направления движения среды. Процесс описывается так называемым волновым уравнением в движущейся среде, получающимся из обычного волнового уравнения и„ -сг (г;Л( + uyv) = 0 при замене координат х' = х-at:

дии + 2адь.и + агдх.х.и - с2 [дх,х.и + д)уи) = 0. (2.1)

Здесь с - скорость распространения волн в покоящейся среде,

О<а <с скорость движения среды, обе величины (с и а) являются постоянными. В полярных координатах уравнение (2.1) имеет более сложный вид, в частности, содержит смешанные производные по времени и пространству.

В. Процесс распространения акустических волн в двухслойной среде описывается волновым уравнением с переменным коэффициентом (отвечающим за скорость распространения волн в среде):

д„и - с2 + + = 0, (2.2)

где с(р) = с, = с, 0<<р<л и с{<р) = с2 =с-а, к<(р<2п.

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

|5„и-Д„« = г, (*,>>) еМ2, />0, (2з)

1ыи="с 5,4=„=г'р

где оператор Дс%а обозначает оператор Лапласа в движущейся или слоистой среде соответственно. Предполагается, что начальные данные и0, и! и правая часть g имеют носитель, находящийся внутри области г<Кг.

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

Работа ведется с дискретным аналогом задачи (2.3), построенным для некоторого, определенного заранее способа дискретизации уравнений (разностная схема) на фиксированных сетках по времени и пространству (Тг и П* соответственно). Обе сетки предполагаются равномерными с шагами г и /г = (/гл,/г?)), при этом ... < г, <ЙГ = г0 <гх <... и Ь^ = 2/т / М, где М - это

количество точек сетки О? по переменной <р. Для дискретизированных операторов Э„ и Ас а вводятся обозначения Оти и , дискретное решение записывается как и'ь, а его значения в точке / = 1р, г = г^ <р = <рт обозначается через ир1т. Таким образом, дискретный аналог дифференциальной задачи (2.3) имеет вид:

р>у*-д:у*Р> о,

[<„,= О, 0, 7 = ...,-1,0,1,...; т = 0,...,М-I, р = 0. (2А)

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

Доказывается, что точные ДПГУ можно сконструировать из решений набора дискретных внешних начально-краевых задач, имеющих вид:

[РФ-Х^У =0, р>0, у>0, /л = 0,...,М-1,

(с1"')" =0, =0, т = 0,...,Му>0, (2.5)

V /),т \ / 7,/н

(С1"')'

V /0 т

т '

Здесь 5Р - это дискретная дельта функция по времени, <5° = 1, 8Р = 0 при р> 0; Ф(п), /? = 0,1,...,М-1 - это некоторый полный дискретный базис на границе г = Лг (базис из синусов и косинусов обычно является удовлетворительным). Стоит отметить, что решения задач (2.5) есть не что иное, как проекции на базис дискретной функции Грина внешней начально-краевой задачи Дирихле. Показывается, что точные ДПГУ определяются соотношением между значениями решения при г = Яг и при г = Яг+кг. Действие ДПГУ заключается в разложении решения на границе г = Яг по базису дальнейшей свертки каждого элемента разложения с соответствующей «базисной» дискретной функцией Грина и обратной сборке решения по базису:

(2.6)

р-' 6*р >1 '

I Лг+Л,

Соотношение (2.6) используется для обновления решения на вспомогательной границе г = Яг + Иг при наличии найденного по разностной схеме для внутренней области (г<Яг) решения при г = Яг. В формуле (2.6) использованы обозначения Р, Г"1 для прямого и обратного разложения по базису ф'"', символ * обозначает дискретную свертку по времени, а матрица ¿> определяется через коэффициенты разложения с'"1 по базису Ф("'

и-I .

на вспомогательной границе г = Яг , т.е. ш^Л Ф(Л). Таким об-

1 А=0

разом

С =

¿0.0

00..

с'-1

¿1М-,

(2.7)

где элементы С"'1' зависят от временного слоя ¡р.

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

\i-n„

Рис 1 Схематичное изображение аппроксимации матрицы G по upociранетву

Предлагается аппроксимировать матрицу G (которая в случае а Ф О является полностью заполненной) урезанной матрицей G с элементами G"'k при п<п1ВС и \rt-к\< пЬат1, при этом все остальные элементы полагаются равными нулю (см. Рис. 1). Параметры nusc и подбираются на практике экспериментально.

Возможность пренебрежения элементами при п > п1ВС основано на предположении гладкости решения, что позволяет считать вклад высоких гармоник по угловой переменной небольшим. Исключение элементов вне ленты jn-/(|<nw предполагает малость этих элементов, наблюдающуюся при умеренной анизотропии (для изотропии эти элементы равны нулю, поскольку матрица G диагональна).

В параграфе 2.4 описывается аппроксимация оператора ДПГУ по времени. основанная на алгоритме, предложенным A. Arnold, М. Ehrhardt, I. Sofronov (2003) для построения дискретных граничных условий для уравнения Шредингера.

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

(2.8)

Здесь а"* - коэффициенты разложения, q" k отвечают за показатели экспонент, L"'k- количество экспонент (квадратные скобки используются для обозначения возведения в степень).

Использование аппроксимации (2.8) позволяет вычислять свертку в ДПГУ (2.6) по рекуррентным формулам

МЧ^Ч'/Г. (/;')-'= 0.

(2.9)

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

дополнительными «накопительными» переменными 1",к. Такая численная локализация нелокальных ДПГУ существенно сокращает объемы вычислений и является ключевой составляющей метода.

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

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

ю" ю': 10"' 10"' 10"

•

в ......"А "^л В "

A "I ...

100 200 300 " 0 100 200 300 400 500

р - number of the time step p - number of the time step

(a) (6)

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

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

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

• решение UDTBC, полученное с использование ДПГУ на г = Rr;

• решение UEXT, полученное на той же сетке что и UDTBC, но в расширенной области, так что отражения от внешней границы не достигают области г < Rr за рассматриваемое время моделирования;

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

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

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

5иотвс (0 = Ц^нг -^вгасЦ^+ытЦ-Совместная ошибка схемы и граничных условий есть

^агрг+ОТБС (0 =

Ошибка аппроксимации схемы определяется как

Оператор ДПГУ не ухудшает точность аппроксимации разностной схемы, если максимальная (по времени) ошибка 311 +отвс{/) близка к ошибке

81]арр1. В частности, если 5иотвс существенно меньше погрешности аппроксимации разностной схемы (например, более чем на порядок), то условие бит„опс « биаррг выполнено.

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

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

Ресурсы, необходимые для вычисления и аппроксимации матрицы <2, велики, однако эти вычисления достаточно провести только один раз (сетка, разностная схема, внешняя граница г = Яг, параметры среды для внешней области фиксированы), а затем использовать для расчетов множества внутренних задач. При этом допускается изменение любых параметров внутренней задачи: источников, начальных данных и даже уравнений и их параметров, при условии, что изменение сосредоточено строго внутри г<Яг.

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

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

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

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

Для демонстрации точности ДПГУ проводятся расчеты для набора параметров птвс и пЬаП(1. Анализ результатов подтверждает, что при достаточно больших значениях п1вг и пыы ДПГУ не ухудшают точность разностной схемы. При этом в рассматриваемых экспериментах для самой подробной сетки (63) точность разностной схемы ~ I О-3, а для самой грубой сетки (С1) точность порядка пары процентов. Точность ДПГУ 81]В1ВС ограничена снизу погрешностью аппроксимации С! по времени. Приводятся численные эксперименты, демонстрирующие, что точность аппроксимации улучшается при увеличении количества экспонент Ь.

(а) _ (б)

Рис 3 Нормы элементов матрицы 6 в зависимости от номера диагонали п-к для сетки 01. (а) а = 0 2, (б) а = 01

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

новная ошибка при аппроксимации 0"л суммами экспонент, не присутствуют в точном решении.

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

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

4

-20

-10

10

20

0

п-к

(б)

Рис. 4 Точность аппроксимации СГ0(.£ = 40, круглые значки) и С'"°(1 = 100, треугольные значки), для каждого элемента б"1 как функция номера диагонали п-к для сетки 03 (а) Ь, -норма по времени (совпадает с нормой во всем частотном диапазоне) (б) Точность после исключения высоких частот

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

В параграфе 3.1 описывается основная идея метода, поясняются его преимущества и недостатки.

Предложенные в предыдущей главе ДПГУ зависят от разностной схемы, что обусловлено наличием дискретной дельта функции по времени в краевых условиях задач (2.5). Для исключения этой зависимости предлагается заменить дискретную дельта функцию 5Т на ее гладкий аналог 5'^). Требуется, чтобы с помощью 51 (/) можно было бы аппроксимировать гладкие функ-

Рис 5. Локальный сплайн 8] (г). Точками помечены значения функции на сетке Т'.

ции с порядком не меньшим чем порядок разностной схемы относительно т. Гладкость ¿¡(f) позволяет перейти от дискретной постановки к дифференциальной и тем самым обеспечить независимость граничных условий от разностной схемы (зависимость от параметров сетки Л, и г сохраняется). Дополнительно G,, в силу гладкости, аппроксимируется меньшим количеством экспонент, что уменьшает расходы на реализацию граничных условий, однако приводит к возрастанию предварительных расходов на вычисление G(.

В параграфе 3.2 описываются детали метода построения ДПГУ со сглаженной функцией Грина. Построение проводится на сетках V и Q*, определенных в главе 2. В качестве <?sr(t) используется локальный сплайн

B.C. Рябенького, построенный по значениям дискретной ST (см. Рис. 5). Такая функция непрерывна вплоть до 3-ей производной, имеет локальный носшель [-2г,2г] и обеспечивает аппроксимацию произвольной достаточно гладкой функции с порядком г'.

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

5„G'M-

G,

-г.(>01

- дс „с; = о,

г и —¿г

I > -2т, г > R.,

-0,

(3.1)

Доказывается, что для гладкой функции /, являющейся суммой не более М гармоник Ф'^ по переменной <р, решение внешней начально-краевой

задачи

б „и - Дь аи =0, / > 0, г > /?,,

«Ц =/м

аппроксимируется с порядком г' на слое /- = /■; дискретной функцией йгН = *Р^/Г,/1 Делается вывод, что использование ДПГУ вида

и'"] =Г

Ц -А,

;>*F[}rV

не нарушает второй порядок аппроксимации

разностных схем.

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

Для нахождения сглаженной функции Грина задачи (3.1) решаются численно по разностной схеме, используемой в дальнейшем и для решения внутренней задачи (центральные разности второго порядка). Для корректного описания гладкой 8[ (/) сетка 63 измельчается во всех направлениях

в с! раз. Элементы сглаженной приведены на Рис. 6 для двух значений с!, (1 = 4 и ¿У = 8. Вычисленные 6"'к являются гладкими функциями с незначительными осцилляциями, вызванными недостаточной точностью решения задач (3.1), при измельчении сетки (увеличении ¿/) амплитуда ос-цилляций уменьшается.

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

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

ю

ю

ю

-DGF -MDGF, 4 -MDGF, 8

«

:

05

1 5

time

(а)

time

(б)

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

га Ю

2 10 Е

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

( 1 Л2 ) д и - дц+—ди—-и =0,Л = 2л-т (3.2)

I г г )

Сравнение проводится на примере одного параметра т (т = 6).

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

Анализ результатов показывает, что для рассматриваемой задачи, несмотря на спектральную скорость сходимости, метод, использующий функции Лагерра и псевдоспектральную дискретизацию по пространству, является более ресурсоемким, чем обычная явная центрально-разностная аппроксимация второго порядка вплоть до точностей КГ6 (см. Рис. 7).

Четвертая глава посвящена аналитическим прозрачным граничным условиям (ПГУ).

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

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

и: г л _ [_ г с. : —*-РО " -®"5РЕСТРА1.

==:\

-1 1-1 г п---г - -1 п ---1-4 -1 т 1Т - - г -,

м 1" Ч"" 1~!" 1" ■¡п [ \ 1 1 :л"гг " ч 1 1 1-г - *

- 1ч » М - - - - Г .- 1 и у * 1 -

Г, | ,1,1 V - - *

-- - Нг ; !

10'

10" 1(Г 1<Г

Рис 7 Количество операций относительно точности решения Сплошная кривая соответствует конечно-разностному подходу, пунктирная - спектральному

Рис 8 Расчетная область

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

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

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

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

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

ди, да„ сг ег Зсг

Р-=- =-+

Ы дг г г дг

ди, Эсг а„ да,, р—- = ——+ -

(4.1)

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

дагг _ диг иг диг

~~эГ" '2 + Л13"Т~>

д( дг г дг ^-А ^+ А ^+ А ^

- /1| I /1,1 Т

5/ 12 дг 11 г 13 дг

5/ ~Аз дг 13 г А" & '

(4.2)

д1 ц\дг дг

Здесь (и, ,и.) и (ст|;.,ст;Л,,<т„,сг.;) - это компоненты вектора скорости и тензора напряжений соответственно, р - плотность, А1,, Ап, Аи, А33, А44 -параметры среды.

Рассматривается процесс распространения волн в бесконечной полуполосе, 0 < г < +со, <2< 2ПИХ, с граничными условиями:

и2 = 0, а,. = О,

на верхней и нижней границах (: = 2тт, г = /Гта1(). Вычислительной областью является ограниченный прямоугольник [О, /?г] х [2т1П,], на правой границе которого г = Лг строятся прозрачные граничные условия, обеспечивающие прохождение волн без отражений.

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

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

д_

дг

= Р(г,з,1,р,Лп,Лв,Л3},А44)

(4.3)

где иг, - это иг, агг компоненты после преобразования Фурье по переменной - и преобразования Лапласа по переменной /, а $ и / - соответствующие двойственные переменные. Матрица

Р =

(х+-х~)

-2

-г'С-1'3 ■ , ,1 2 -1--/С,/М/ -Я,,/ -Я2/05

V 5

2 0 1;

где скалярные функции х± определены как

===== К\{г-\а1г + /)р$1 + № + ПР^? + V

К, ^а!2 + /V + + црь1!"- + V

Здесь К{ - функция Макдональда первого порядка, а параметры <у0, а2, , с,, с.,, уЗ, 77, зависят только от параметров среды:

а„ = —-а2

2АпА4А

Д ] + Д4 ^ ( Дц + Дф

2ЛПД

44

— ~ ; > с-1 — 2Д4а0,

2Д,Д,4 2Д,Д44 Д, ДЦ^44

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

1 „ л, \ 1

Г 5

(4.4)

матрицы £>0, Qí, ()2 зависят только от параметров среды, матрицы К (5) зависят дополнительно как от номера Фурье гармоники, так и от положения границы Элементы матрицы К (я) стремятся к 0 при .у —>+со. После обратных преобразований Фурье и Лапласа условий (4.3) с использованием представления (4.4) ПГУ в физических переменных принимают вид:

а

д_

дг

+ ~во

г

+а Ч

(Т.,

(4.5)

4т = 0.

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

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

В проведенных численных экспериментах граничные условия (4.5) используются совместно со схемой .1. Ушеих (1986) на сдвинутых сетках для системы (4.1), (4.2). На границе вычислительной области г = Кт компоненты иг и ап обновляются с использованием дискретизированных уравнений (4.5), в то время как для обновления остальных компонент используются дискретизированные уравнения основной системы (4.1), (4.2).

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

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

(а) (б)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2. Sofronov I.L, Podgornova О. 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. Подгорнова OB Построение оператора дискретных неотражающих граничных условий для моделирования волн в движущейся среде // Математическое моделирование - 2007. - Т. 19, № 8. - С. 75-82.

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

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

6. Podgornova О. V., Sofronov I. L. Toward efficient numerical generation of low-reflecting boundary conditions for anisotropic media // Abstracts. International seminar «Days on diffraction» - Saint Petersburg, 2005. -P. 73-74.

7. Подгорнова O.B., Софронов И.Л. Эффективное вычисление и аппроксимация функций Грина для оператора слабо-отражающих гра-

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

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

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

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

Автореферат

Заказ № 92/09/08 Подписано в печать 11 09 2008 Тираж 80 экз Уел пл 1,25

/: - "о ООО "Цифровичок", тел (495) 797-75-76, (495) 778-22-20 www.cfr.ru; е-тш1:т/о@с/г.ги

ВВЕДЕНИЕ

ГЛАВА 1. МЕТОДЫ ПОСТРОЕНИЯ НЕОТРАЖАЮЩИХ ГРАНИЧНЫХ

УСЛОВИЙ

1.1 Понятие неотражающих граничных условий

1.2 Обзор методов

1.2.1 Локальные условия.

1.2.2 Идеально поглощающий слой (PML).

1.2.3 Нелокальные условия.

1.3 Пример построения прозрачных граничных условий

ГЛАВА 2. ДИСКРЕТНЫЕ ПРОЗРАЧНЫЕ ГРАНИЧНЫЕ

УСЛОВИЯ (ДПП/)

2.1 Модельные задачи

2.2 Точные ДПП/

2.3 Пространственная аппроксимация ДПГУ

2.4 Аппроксимация ДПГУ по времени

2.4.1 Аппроксимация суммой экспонент.

2.4.2 Алгоритм построения аппроксимации.

2.4.3 Численный пример.

2.5 Точность ДПГУ: типичное поведение ошибки

2.6 Вычислительные затраты на ДПГУ

2.7 Численное исследование: волновое уравнение в движущейся среде

2.7.1 Дискретизация задачи.

2.7.2 Параметры численного эксперимента.

2.7.3 Результаты применения ДПГУ.

2.7.4 Расчеты на большие времена.

2.7.5 Апостериорная оценка точности.

2.8 Численное исследование: волновое уравнение в слоистой среде

2.9 Обсуждение результатов

ГЛАВА 3. ДПП/ СО СГЛАЖЕННЫМ ЯДРОМ

3.1 Описание задачи

3.2 Построение ДПГУ

3.3 Численное исследование

3.3.1 Параметры численного эксперимента.

3.3.2 Ядро ДПГУ.

3.3.3 Результаты применения ДПГУ.

3.3.4 Выводы.

3.4 Сравнение эффективности конечно-разностного и спектрального методов для вычисления ядра ДПГУ

3.4.1 Постановка задачи.

3.4.2 Спектральная дискретизация.

3.4.3 Свойства разложений по функциям Лагерра.

3.4.4 Алгоритм решения в терминах коэффициентов ряда Лагерра.

3.4.5 Сравнение эффективности.

3.4.6 Выводы.

ГЛАВА 4. ДВА ПРИМЕРА АНАЛИТИЧЕСКИХ ПРОЗРАЧНЫХ

ГРАНИЧНЫХ УСЛОВИЙ

4.1 Аналитические ПГУ для одной задачи аэроакустики

4.1.1 Постановка задачи.!.

4.1.2 Неотражающие граничные условия.

4.1.3 Численные результаты.

4.2 Аналитические ПГУ для линейной упругости в анизотропной среде

4.2.1 Постановка задачи.

4.2.2 Построение ПГУ.

4.2.3 Асимптотическое разложение граничного оператора.

4.2.4 Свойства ПГУ.

4.2.5 Частный случай ПГУ для изотропной среды.

4.2.6 Аппроксимация ПГУ.

4.2.7 Дискретизация ПГУ.

4.2.8 Численные эксперименты.

Актуальность темы

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Теоретическая и практическая ценность

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

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

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

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

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

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

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

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

Публикации

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

В работах с соавторами лично соискателем выполнено следующее: [18,60,10] - окончательно сформулирован метод, разработан и реализован алгоритм построения дискретных прозрачных граничных условий, проведены численные эксперименты на двух модельных задачах; [10,56] - предложен и реализован способ конструкции дискретных прозрачных граничных условий, непривязанный к численному методу решения основной задачи; [19] - скорректирован алгоритм построения аппроксимации дискретных граничных условий, получены результаты лучшей точности; [55] — проведено сравнение спектрального и конечно-разностного подходов для нахождения дискретной функции Грина внешней задачи; [2] - реализован оператор прозрачных граничных условий для линеаризованной системы уравнений Эйлера, проведены численные эксперименты.

Апробация

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

• международной конференции «Workshop on nonlinear approximations in numerical analysis» (Москва, 2003),

• XV и XVI Всероссийской конференции «Теоретические основы и конструирование численных алгоритмов и решение задач математической физики с приложением к многопроцессорным системам», посвященная памяти К.И. Бабенко (Абрау-Дюрсо, 2004, 2006),

• На международном семинаре «Days on diffraction» (Санкт-Петербург, 2005),

• на Всероссийской научно-практической конференции «Вычислительный эксперимент в аэроакустике» (Светлогорск, 2006),

• на Всероссийской конференции по вычислительной математике «КВМ-2007» (Новосибирск, 2007),

• на международной конференции «Matrix methods and operator equations» (Москва, 2007).

Структура и объём диссертации

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

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

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

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

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

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

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

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

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

Заключение

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

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

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

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

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

Мы численно исследуем ДПГУ на двух модельных двумерных задачах: волнового уравнения в движущейся среде и волнового уравнения в слоистой среде. Первая задача, которой посвящен параграф 2.7, представляет собой пример анизотропной среды, вторая, описанная в параграфе 2.8, - неоднородной среды. Мы демонстрируем высокую точность построенных граничных условий (~1(Г3) и их устойчивость на больших временах моделирования при помощи представительного набора численных экспериментов на вложенных сетках. Мы показываем, что основным параметром, определяющим точность граничных условий, является количество экспонент, используемых для аппроксимации, что для нашего алгоритма эквивалентно временному интервалу, на котором производится интерполяция.

В дополнении к основному подходу построения ДПГУ в главе 3 предлагается метод, основанный на сглаженной функции Грина. Отличием от подхода предыдущей главы является изменение внешних начально-краевых задач, решения которых определяют ядро ДПГУ (функцию Грина). В подходе главы 3 дискретные внешние задачи заменяются дифференциальными, со специальными, гладкими краевыми условиями на внутренней границе. Для описания зависимости краевых условий по времени мы используем аппроксимацию дискретной дельта-функции, полученную с помощью локального сплайна. В параграфе 3.2 мы доказываем, что ДПГУ со сглаженным ядром сохраняют второй порядок аппроксимации разностных схем. Достоинством таких ДПГУ является независимость численных методов для решения внутренних и внешних задач, а к недостаткам относится большее количество вычислительных ресурсов, требуемых для нахождения функции Грина. В параграфе 3.3 мы приводим численные результаты, демонстрирующие высокую точность подхода (также вплоть до 1(Г3).

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

Глава 4 посвящена аналитическим прозрачным граничным условиям для двух задач: системы уравнений Эйлера и уравнений линейной упругости в анизотропной среде. В параграфе 4.1 мы применяем прозрачные граничные условия для линеаризованной системы уравнений Эйлера [28] в тестовой аэроакустической задаче являющуюся моделью звукопоглощающей конструкции. Мы сравниваем точность ПГУ, нелокальных, построенных для линеаризованных уравнений и точность локальных граничных условий для исходной нелинейной системы [7]. Проведенные эксперименты демонстрируют лучшую точность ПГУ.

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

1. Абрамовиц М., Стиган И. (Ред.) Справочник по специальным функциям с формулами, графиками и таблицами. Пер. с англ. М.: Наука, 1979. -832с.

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

3. Бабенко К.И. и др. Пространственное обтекание гладких тел идеальным газом. М.:Наука, 1964. - 508 с.

4. Бейкер Дж., Грейвс-Моррис П. Аппроксимации Паде. Основы теории. Обобщения и приложения. М.: Мир, 1986. - 502 с.

5. Владимиров B.C. Уравнения математической физики. М.:Наука, 1976. -527 с.

6. Дёч Г. Руководство к практическому применению преобразования Лапласа и Z-преобразования. М.: Наука, 1971. - 288 с.

7. Дородницын JI.B. Неотражающие граничные условия для систем уравнений газовой динамики // Ж. вычисл. матем. и матем. физ. 2002. - т.42, №4. - С. 522-549.

8. Лисица В.В. Спектрально-согласованные сетки для моделирования волновых процессов: Диссертация на соискание ученой степени кандидата физико-математических наук. СО РАН, Институт Нефтегазовой Геологии и Геофизики им. А.А. Трофимука, 2007. - 90 с.

9. Лурье А.И. Теория упругости. М.:Наука, 1970. - 940 с.

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

11. Рябенький B.C. Введение в вычислительную математику. -М.:Физматлит, 2000. 294 с.

12. Рябенький B.C. Метод разностных потенциалов и его приложения. -М.: Физматлит, 2002. 420 с.

13. Рябенький B.C. Точный перенос разностных краевых условий // Функ. анализ и его приложения. 1990. - Т.24, В.З. - С. 90-91.

14. Софронов И.Л. Условия полной прозрачности на сфере для трехмерного волнового уравнения // Доклады РАН. 1992. - Т. 326, №6. - С. 453-457.

15. Софронов И.Л. Условия полной прозрачности для волнового уравнения: Препринт / ИПМ им. М.В. Келдыша РАН. М., 1993. - №76. - 25 с.

16. Софронов И.Л. Точные искусственные граничные условия для некоторых задач аэродинамики и дифракции: Диссертация на соискание ученой степени доктора физико-математических наук. — Институт Прикладной Математики им. М.В. Келдыша РАН, 1999. -215 с.

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

18. Суетин П.К. Классические ортогональные многочлены. — М.'Физматлит, 2005.-480с.

19. Abarbanel S., Gottlieb D. A mathematical analysis of the PML method // J. Comput. Phys. -1997. 134, P. 357-363.

20. Abarbanel S., Gottlieb D. On the construction and analysis of the absorbing layers in СЕМ // in 13th Annual Review of Progress in Applied Computational Electromagnetics. 1997. - P. 876-883.

21. Alpert В., Greengard L., Hagstrom T. Rapid evaluation of nonreflecting boundary kernels for time-domain wave propagation // SIAM J. Num. Anal. -2000.-37, n4.-P. 1138-1164.

22. Alpert В., Greengard L., Hagstrom T. Nonreflecting Boundary Conditions for the Time-Dependent Wave Equation // Journal of Computational Physics.2002.- 180.-P. 270-296.

23. Arnold A., Ehrhardt M., Sofronov I. Discrete transparent boundary conditions for the Schrodinger equation: Fast calculation, approximation, and stability // Comm. Math. Sci. 2003. - 1, No.3. - P. 501-556.

24. Asvadurov S., Druskin V., Guddati M.N., Knizhnerman L. On optimal finite-difference approximation of PML // SIAM Journal of Numerical Analysis.2003.-41.-P. 287-305.

25. Asvadurov S., Druskin V., Knizhnerman L. Application of the difference Gaussian rules to solution of hyperbolic problems II. Global expansion // Journal of Computational Physics. 2002. - n. 175 - P. 24-29.

26. Ballmann J., Britten G., Sofronov I. Time-Accurate Inlet and Outlet Conditions for Unsteady Transonic Channel Flow // AIAA Journal. 2002. - Vol. 40, No. 9.-P. 1745-1754.

27. Bayliss A., Turhel E. Radiation boundary conditions for wave-like equations // Comm. Pure Appl. Math. 1980. - 33. - P. 707-725.

28. Becache E., Fauqueux S., Joly P. Stability of Perfectly Matched Layers, Group Velocities and Anisotropic Waves // J. Comput. Phys. 2003. - 188 (2). - P. 399-433.

29. Becache E., Petropoulos P.G., Gedney S.D. On the Long-Time Behavior of Unsplit Perfectly Matched Layers I I IEEE Trans. Antennas Propagat. 2004. -52(5).-P. 1335-1342.

30. Berenger J.-P. A perfectly matched layer for the absorption of electromagnetic waves // J. Comput. Phys. 1994. - 114. - P. 185-200.

31. Boyd John P. Chebyshev and Fourier Spectral Methods. New York: DOVER Publications, Inc., 2001.-691 p.

32. Collino F., Tsogka C. Application of the PML absorbing layer model to the linear elasodynamic problem in anisotropic heterogeneous media: INRIA, Rapport de recherche, 1998. n° 3471. - 29 p.

33. Chung Young-Seek, Sarkar Tapan K., Jung Baek Ho, Salazar-Palma Magda-lena An Unconditionally Stable Scheme for the Finite-Difference Time-Domain Method // IEEE Transactions on Microwave Theory and Techniques. 2003. - Vol. 51, No. 3. - P. 697-704.

34. Dedner A., Kroner D., Sofronov I.L., Wesenberg M. Transparent Boundary Conditions for MHD Simulations in Stratified Atmospheres // J. Comput. Phys. 2001. - 171. - P. 448-478.

35. Diaz J., Joly P. An analysis of higher order boundary conditions for the wave equation // SIAM J.Math. 2005. - 65. - P. 1547-1575.

36. Diaz J., Joly P. A time domain analysis of PML models in acoustics // Comput. Meth. Appl. Mech. Eng. 2006. - 195. - P. 3820-3853.

37. Engquist В., Majda A. Absorbing boundary conditions for the numerical simulation of waves // Math. Сотр. 1977. - 31. - P. 629-651.

38. Givoli D. Nonreflecting boundary conditions // J. Comput. Phys. 1991.-94, l.-P. 1-29.

39. Givoli D., Cohen D. Non-reflecting boundary conditions based on Kirchhoff-type formulae // J. Comput. Phys. 1995. - 177. - P. 102-113.

40. Grote M.J., Keller J.B. Exact nonreflecting boundary conditions for the time dependent wave equation // SIAM J.Appl.Math. 1995. - 55. - P. 280-297.

41. Grote M.J., Keller J.B. Nonreflecting boundary conditions for Maxwell's equations // Journal of Computational Physics. 1998. - v. 139, n.2. - P. 327342.

42. Grote M.J., Keller J.B. Exact nonreflecting boundary conditions for elastic wave //. SIAM Journal of Applied Mathematics. 2000. - v.60, n.3. - P. 803819.

43. Grote M.J., Kirsch C. Nonreflecting boundary condition for time dependent multiple scattering // J. Comput. Phys. 2007. - Vol. 221, l.-P. 41-62.

44. Hagstrom T. Radiation boundary conditions for the numerical simulation of waves // ActaNumer., 1999. 8, 47. - P. 47-106.

45. Hagstrom Т., Hariharan SJ., Thompson D. High-Order Radiation Boundary Conditions for the Convective Wave Equation in Exterior Domains // SIAM J. Sci. Comput., 2003.-Vol. 25, No. 3.-P. 1088-1101.

46. Hagstrom Т., Mar-Or A., Givoli D. High-order local absorbing conditions for the wave equation: Extensions and improvements // J. Comput. Phys., 2008. -227.-P. 3322-3357.

47. Higdon R.L. Numerical absorbing boundary conditions for the wave equation // Math. Comput., 1987. 49. - P. 65-90.

48. Wolfram MathWorld, a free resource from Wolfram Research built with Mathematica technology, <http://mathworld.wolfram.com> (2008, September).

49. Maplesofit, a division of Waterloo Maple Inc. <http://www.maplesoft.com/>, (2008, September).

50. Konyitkh G.V., Mikhailenko B.G., Mikhailov A.A. Application of the Integral Laguerre Transforms for Forward Seismic Modeling // Journal of Computational Acoustic, 2001. Vol. 9, No. 4. - P. 1523-1541.

51. Lyrintzis A.S. Review: The use of Kirchhoffs method in computational aeroacoustics // Journal of Fluids Engineering, 1994. 116 №12. — P. 665676.

52. Mikhailenko B.G., Mikhailov A.A., Reshetova G.V. Numerical viscoelastic modeling by the spectral Laguerre method // Geophysical Prospecting, 2003. -51.-P. 37-48.

53. 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 and operator equations» -Moscow, 2007. C. 66-67.

54. Podgornova O.V., Sofronov I. L. Toward efficient numerical generation of low-reflecting boundary conditions for anisotropic media // Abstracts. International seminar «Days on diffraction» Saint Petersburg, 2005. — P. 73-74.

55. Ryaben'kii V.S., Tsynkov S.V., Turchaninov V.I. Global Discrete Artificial Boundary Conditions for Time-Dependent Wave Propagation // J.Comput. Phys., 2001. 174. - P. 712-758.

56. Sofronov I.L. Artificial boundary conditions of absolute transparency for two-and three-dimensional external time-dependent scattering problems // Euro. J.Appl.Math., 1998. -V.9, No.6. P. 561-588.

57. Sofronov I.L. Non-reflecting inflow and outplow in wind tunnel for transonic time-accurate simulations // J. Math. Anal. Appl. 1998. - 221. - P. 82-115.

58. Sofronov I.L., Podgornova О. V. A spectral approach for generating non-local boundary conditions for external wave problems in anisotropic media // Journal of Scientific Computting. 2006. - V. 27, No 3. - P. 419-430.

59. Sofronov I., Zaitsev N. Non-reflecting boundary conditions for anisotropic media // Book of Abstracts. 77 Annual Meeting of the GAMM. -TU Berlin, 2006. P. 362.

60. Taflove A. Computational Electrodynamics: The Finite-Difference Time-Domain Method. Boston: Artech House, 1995. 599 p.

61. Тат C.K. W., Auriault L., Cambuli F. Perfectly matched layer as an absorbing boundary condition for the linearized Euler equation in open and ducted domains // J. Comput. Phys. 1998. - 144(1). - P. 213-234.

62. Ting L., Miksis M.J. Exact boundary conditions for scattering problems // J. Acoust. Soc. Amer. 1986. - 80. - P. 1825-1827.

63. Thomsen L. Weak elastic anisotropy. // Geophysics. 1986. - vol. 51(10). -P.1954-1966.

64. Trefethen L. Spectral methods in matlab. Philadelphia: SIAM, 2000. 160 p.

65. Trefethen L.N., Halpern L. Well-Posedness of One-Way Wave Equation and Absorbing Boundary Conditions // Mathematics of Computation. 1986. -Vol. 47, No. 176.-P. 421-435.

66. Tsvankin I. Seismic signatures and analysis of reflection data in anisotropic media. Elsevier Science Ltd., 2001. 454 p.

67. Tsynkov S. V. Numerical solution of problems on unbounded domains. A review. // Appl. Numer. Math. 1998. - 27(4). - P. 465-532.

68. Virieux J. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. // Geophysics. 1986. - v.51(4). - P. 889-901.