Численное моделирование волн в неограниченных областях

[Rangok]

Munin в сообщении #548132 писал(а):

А вообще, какая там величина этих остаточных ошибок, по отношению к исходной волне?

Если нагрузку прикладывать к внешней границе области, то все нормально: ошибки начинают уменьшаться при увеличении густоты сетки и ширины pml-слоев. В принципе, не особо напрягаясь можно уменьшить величину накладываемых ошибок до 0,5% от максимальной амплитуды волны.

Munin в сообщении #548132 писал(а):

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

Я если честно пока сам не знаю, какая там погрешность меня может устраивать, а заказчиков у меня нет никаких.

Munin в сообщении #548132 писал(а):

а во-вторых, может быть, она - как раз расчётная для заданных параметров pml-слоя?

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

Munin в сообщении #548132 писал(а):

Действительно, странно, особенно то, что "иголки" появляются ещё до того, как волна успела фактически дойти до pml-слоя. Боюсь, ответ в каких-то нюансах численной схемы, как основной, так и для реализации pml-слоя. Что-то на что-то "незаконно" влияет. Подобные "иголки" типичны для несходящихся численных методов, но ничего конкретнее сказать не могу.

Мне кажется, я понял природу появления этих ”иголок”.
Рассмотрим решение в обычной постановке.
Когда мы прикладываем нагрузку с внешней стороны, то на первых временных шагах у нагружаемой границы формируются пики ошибок. По мере прохождения времени волна проходит вправо вместе с этими маленькими ошибками, которые затем сглаживаются и формируют маленькую волну, которая находится перед основной. Причем эта маленькая волна так и остается в решении (http://savepic.su/1578455.jpg). Если делать сетку гуще, то ошибки уменьшаются. По мере продвижения волны вправо градиент напряжений в элементах, прилегающей к нагружаемой плоскости уменьшается, так как соседние к ним элементы набирают все большее и большее напряжение. В результате меньшего градиента напряжений больше никаких пиков не возникает, и решение становится устойчивым (http://savepic.su/1624534.jpg).
Но, ежели мы прикладываем нагрузку внутри области вдоль одной прямой линии узлов, то в элементах, которые касаются линии нагруженных узлов своей правой стороной, возникает растяжение, а в элементах, которые касаются линии своей левой стороной – сжатие. В результате чего происходит образование двух волн: одна вправо (волна сжатия), другая влево (волна растяжения) . В итоге получается, что на линии нагруженных узлов, напряжения равны нулю, а на соседних линиях напряжения равны друг другу, но противоположны по знаку. И по мере возрастания нагрузки картина не меняется, а лишь усугубляется. На линии нагрузки напряжения так и будут оставаться нулевыми, а на соседних линиях лишь увеличиваться. Все это приводит к образованию гигантского градиента напряжений в элементах, граничащих с линией нагружения, и как результат – приводит к огромному острому пику вдоль линии нагрузки (http://savepic.su/1636822.jpg). Причем, сгущение сетки никак не может улучшить данное стечение обстоятельств (http://savepic.su/1634774.jpg, http://savepic.su/1628630.jpg). Когда же нагрузка полностью снимается градиент напряжений уменьшается и ошибок больше не возникает, а предыдущие ошибки сглаживаются.
Если теперь рассмотреть решение в постановке pml с нагруженной внешней стороной, то мы увидим схожую с обычной постановкой картину (http://savepic.su/1630678.jpg). Однако теперь ошибки формируются не виде одного пика, а сразу нескольких (http://savepic.su/1623510.jpg). Следует отметить, что по мере прохождения времени ошибки не сглаживаются в маленькую волну, как в обычном решении, а растягиваются в виде маленьких иголок на все решение. Причем по абсолютному значению, в конце концов, ошибки меньше чем в обычном решении, хоть они размазаны по всей области.
Теперь если приложить нагрузку внутри области происходит то, что я уже показывал ранее на картинках. Ошибки формируются при нагружении, а затем размазываются на все решение.
Отсюда следует вывод –генерировать нагрузку внутри области решения с помощью лишь одной линии нагружения нельзя – потому что это ведет к слишком большим ошибкам.
Но вот из-за чего решение в постановке pml так себя ведет пока еще вопрос открытый.

Munin в сообщении #548132 писал(а):

В какой степени солвер вам подконтролен? Вы его целиком написали сами, используете написанный коллегой, используете готовый продукт?

Солвер мне подконтролен полностью, кроме решения СЛАУ (здесь я использую солвер Matlab'a). Напсал я его целиком, даже программу решения СЛАУ писал, только использую все равно Matlab'овскую.

Munin в сообщении #548132 писал(а):

Отдельно тот же вопрос по генерации сетки. И чуть поподробнее, как вы pml-слой задаёте?

Программу генерации области решения я написал полностью сам. Однако, писал я ее, когда только начинал изучать программирование (полтора года назад). Писал по средствам метода проб и ошибок, мало что понимая в каких-либо парадигмах программирования и организации программ. Поэтому, она получилась очень плохая – работает медленно, возможности ограничены, дает сбои и т.д. и т.п. Поэтому, когда передо мной встал вопрос о решении задач с использованием pml-слоев, для которых необходимо создавать дополнительные зоны решения я решил не использовать для этих целей свою программу, а использовать готовый генератор сеток Matlab'a. С помощью него я создаю дополнительные зоны pml, на тех границах, в которых мне это необходимо. В принципе, можно целиком использовать генератор Matlab'a для создания всей сетки, но я по старинке все равно использую свою программу, так как в дальнейшем планирую все равно заново написать целиком собственный генератор для своей задачи, включая и генерацию pml-слоев.

Munin в сообщении #548132 писал(а):

Это тяжко. Но может быть, он хотя бы порекомендует своего коллегу, к которому вы могли бы обратиться?

Да не, никто со мной просто так возится, не будет, все-таки вопросы не самые простые возникают, во все надо целиком вникать, а это никому не надо.

-- Чт мар 15, 2012 09:04:42 --

Suvorov.as в сообщении #548463 писал(а):

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

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

Suvorov.as в сообщении #548463 писал(а):

Далее:
Посмотреть как будут работать PML слои в зависимости от их толщины.

Тоже делал. Чем слой толще - тем меньшие ошибки накладываются.

[Suvorov.as]

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

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

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

[INGELRII]

По поглощающим граничным условиям сказать ничего не могу, незнаком. По сути же:

1) Вы исследовали свою схему на сходимость? Такое подозрение, что она сходима только условно. Острые пики, не исчезающие при измельчении сетки - это ж классический индикатор условной сходимости.

2) Усреднять решение не пробовали? Ну да, не совсем это красиво, а что делать. Зато пики будут побеждены. (Да, искусственная вязкость тоже рулит)

[Rangok]

Suvorov.as в сообщении #548479 писал(а):

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

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

Suvorov.as в сообщении #548479 писал(а):

Напрмер в уравнения можно вводить искуственную вязкость пропорциональную второй производной от функции по пространственной координате. Это сгладит выбросы.

Интересно, а можно про это поподробнее где-нибудь почитать? Не изменит ли такая процедура физику процесса, не отдалит ли ее от реальности?

Suvorov.as в сообщении #548479 писал(а):

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

Я именно так и делаю, посмотрите на рисунки. В "нормальных" решениях профиль волны полностью повсторяет профиль нагрузки.

-- Чт мар 15, 2012 17:09:12 --

INGELRII в сообщении #548535 писал(а):

Вы исследовали свою схему на сходимость? Такое подозрение, что она сходима только условно. Острые пики, не исчезающие при измельчении сетки - это ж классический индикатор условной сходимости.

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

INGELRII в сообщении #548535 писал(а):

Усреднять решение не пробовали? Ну да, не совсем это красиво, а что делать.

Нет не пробовал, а как это делать?

[Suvorov.as]

Немного по вашей проблеме есть здесь

Расчет бистатической силы цели сложных многорезонансных оболочек методом конечных элементов

АКУСТИЧЕСКИЙ ЖУРНАЛ, 2011, том 57, № 5, с. 709–716

[Rangok]

Бесплатно где-нибудь можно скачать?

[INGELRII]

Rangok в сообщении #548592 писал(а):

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

Пардон муа, за невнимательность. Сходимость - это хорошо.

Усреднение - это, скажем, получили Вы на каком-то шаге решение. Затем в каждой точке берем, например, среднее арифметическое значений по всем соседним узлам - и говорим, что это и будет "уточненное" решение на данном шаге. Далее применяем разностные схемы непосредственно к этому "уточненному" решению. То, что получилось, тоже усредняем. Ну и так на каждом шаге.

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

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

[Suvorov.as]

INGELRII в сообщении #548614 писал(а):

Внести в них маленькую вязкость - вполне допустимо.

Для такого класса гиперболических задач это не просто допустимо а необходимо.

Прочитайте: Проуч. Вычислительная гидродинамика.
глава про ударные волны

[Rangok]

Спасибо, что разъяснили про усреднение, надо будет попробовать ради интереса))
Можете подробнее рассказать про искусственную вязкость, как конкретно ее делать? Ну или может литературу какую посоветуете.

-- Чт мар 15, 2012 18:35:38 --

Suvorov.as в сообщении #548619 писал(а):

Прочитайте: Проуч. Вычислительная гидродинамика.
глава про ударные волны

Спасибо!

[Munin]

Rangok в сообщении #548474 писал(а):

Я если честно пока сам не знаю, какая там погрешность меня может устраивать, а заказчиков у меня нет никаких.

Для кого-то же вы решаете эту задачу, хотя бы для научрука, или по крайней мере откуда-то взяли её. Вообще, погрешность на уровне 1% - это хорошо. Иногда "очень хорошо, и лучшего желать нельзя". Иногда кому-то хочется 0,1%, но достигаются они с большими затратами. Иногда и 10% бывает приемлемо.

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

Rangok в сообщении #548474 писал(а):

Да не, никто со мной просто так возится, не будет, все-таки вопросы не самые простые возникают, во все надо целиком вникать, а это никому не надо.

Да, вопрос "за что возиться" обычно тоже решает научрук, а если вы даже ему не нужны... :-( Ничего, кроме соавторства в будущих статьях, как я понимаю, вы предложить пока не можете. Впрочем, вопросы ваши пока сравнительно простые, и помочь вам могут даже просто так.

Rangok в сообщении #548592 писал(а):

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

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

Rangok в сообщении #548474 писал(а):

Солвер мне подконтролен полностью, кроме решения СЛАУ (здесь я использую солвер Matlab'a). Напсал я его целиком, даже программу решения СЛАУ писал, только использую все равно Matlab'овскую.

Программу генерации области решения я написал полностью сам. Однако, писал я ее, когда только начинал изучать программирование (полтора года назад). Писал по средствам метода проб и ошибок, мало что понимая в каких-либо парадигмах программирования и организации программ. Поэтому, она получилась очень плохая – работает медленно, возможности ограничены, дает сбои и т.д. и т.п. Поэтому, когда передо мной встал вопрос о решении задач с использованием pml-слоев, для которых необходимо создавать дополнительные зоны решения я решил не использовать для этих целей свою программу, а использовать готовый генератор сеток Matlab'a. С помощью него я создаю дополнительные зоны pml, на тех границах, в которых мне это необходимо. В принципе, можно целиком использовать генератор Matlab'a для создания всей сетки, но я по старинке все равно использую свою программу, так как в дальнейшем планирую все равно заново написать целиком собственный генератор для своей задачи, включая и генерацию pml-слоев.

Хм. У меня возник вопрос, а как себя ведёт Matlab-овский решатель СЛАУ в тот момент, когда решение везде нулевое или почти нулевое. Вообще, численные методы СЛАУ - отдельная наука по отношению к науке численных методов ДУЧП, и там есть свои источники погрешностей и ошибок. Тут может быть что угодно, например, решатель СЛАУ переключается с одного метода решения на другой (поскольку шибко интеллектуальный).

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

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

[Rangok]

Munin в сообщении #548690 писал(а):

Хм. У меня возник вопрос, а как себя ведёт Matlab-овский решатель СЛАУ в тот момент, когда решение везде нулевое или почти нулевое. Вообще, численные методы СЛАУ - отдельная наука по отношению к науке численных методов ДУЧП, и там есть свои источники погрешностей и ошибок. Тут может быть что угодно, например, решатель СЛАУ переключается с одного метода решения на другой (поскольку шибко интеллектуальный).

Я сомневаюсь, что тут могут крыться какие-то подводные камни, по той простой причине, что решение СЛАУ очень легко проверить на правильность. Нужно просто помножив вектор решения на матрицу иситемы и вычев из этого правую часть. Я проверял, все решается отлично, ошибки на несколько порядков меньше искомых величин.

Munin в сообщении #548690 писал(а):

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

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

Munin в сообщении #548690 писал(а):

Ещё мысль, что pml-слой чисто физически устроен так, что в нём возбуждения умеют распространяться быстрее звука. Это надо его уравнения посмотреть внимательней.

Уравненеия, которые находятся в pml-слоях - это те же уравнения, что и в обычной области решения, только в которых координата, по которой происходит поглощение волн, заменена следующим образом: $x \rightarrow \tilde{x}=\int_{0}^{x}\lambda (x)dx$, а частные производные заменяются на: $\frac{\partial }{\partial x}\rightarrow \frac{\partial }{\partial \tilde{x}}=\frac{1}{\lambda (x)}\frac{\partial }{\partial x} $, где$\lambda (x)=1+f^{e}(x)-i \frac{c_{s}f^{p}(x)}{\omega }$ - "растягивающая" координату комплексная функция. Более подробно об этом можно почитать здесь http://scholarspace.manoa.hawaii.edu/bitstream/handle/10125/20488/M.S.Q111.H3_4091_r.pdf?sequence=1

Munin в сообщении #548690 писал(а):

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

Я столько времени убил на эти pml-слои, а теперь на помойку выкидывать? Да и не думаю, что простейшая сумма Неймана и Дирихле будет лучше.

[Munin]

Rangok в сообщении #548936 писал(а):

Я столько времени убил на эти pml-слои, а теперь на помойку выкидывать? Да и не думаю, что простейшая сумма Неймана и Дирихле будет лучше.

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

[Suvorov.as]

Как там продвигается решение задачи?

[Rangok]

Нормально, вроде бы. Сделал объемный источник волн, в котором нагрузка плавно приложена, и все стало отлично. Волна довольно ровная, никаких "иголок" больше нет, а остаточные ошибки, накладываемые на все решение около 1-2% от максимальной амплитуды волны. Сейчас пытаюсь найти литературу по искусственной вязкости, применительно к динмаическим задачам теории упругости.

[Munin]

Здорово! Порадуйте картинками :-)