Квазиньютоновский метод для неявного связывания сопряженных задач FSI Текст научной статьи по специальности «Математика»

УДК 519.63+532+531.3

КВАЗИНЬЮТОНОВСКИЙ МЕТОД ДЛЯ НЕЯВНОГО СВЯЗЫВАНИЯ СОПРЯЖЕННЫХ ЗАДАЧ Е81

КОПЫСОВ СП., КУЗЬМИН ИМ., НЕДОЖОГИН Н.С., НОВИКОВ А.К., РЫЧКОВ ВН., ТОНКОВ Л.Е.

Институт механики Уральского отделения РАН, 426067, Ижевск, ул. Т. Барамзиной, 34

АННОТАЦИЯ. Рассматривается разделенный подход при решении жесткосвязанных задач взаимодействия гидро/газодинамических потоков и деформируемых тел на нестыкующихся сетках. Предложен алгоритм для решения нелинейной интерфейсной системы уравнений на границе взаимодействия. Представлены результаты моделирования в рамках жесткосвязанного разделенного подхода некоторых тестовых задач.

КЛЮЧЕВЫЕ СЛОВА: связанные задачи, нелинейные системы уравнений, метод Бройдена, ускорение сходимости, моделирование взаимодействие деформируемых тел с потоками.

ВВЕДЕНИЕ

Выделяют два основных подхода связывания задач взаимодействия потоков жидкости (газа) и деформируемых тел (FSI - Fluid and Structure Interaction): монолитный и последовательный [1].

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

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

нагрузки на эти границы получаются разными со стороны конструкции и газа. Алгоритмы, обеспечивающие сопряжение между связанными задачами, например, гидродинамической задачей (CFD) и задачей деформирования конструкции (CSD), подразделяют по степени сопряжения на слабосвязанные/явные (wFSI) и жесткосвязанные/неявные (sFSI). В первом случае обмен данными между задачами выполняется явно, нагрузка для задачи деформирования на шаге t зависит от давления жидкости на шаге t + Dt. То есть, в слабосвязанном подходе, в пределах одного временного шага, обмен граничными значениями происходит единожды, и нет никаких приближений выполняемых отдельными задачами (рис. 1, а). Таким образом, реализуется так называемое «шахматное связывание», при котором программы обоих приложений не обязательно должны одновременно выполняться на вычислительном кластере.

В жесткосвязанном подходе выполняется несколько обменов на каждом временном

шаге и, тем самым, реализуется неявное сопряжение, обеспечивающее более сильную связь

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

выполняют решение нелинейной интерфейсной системы уравнений (рис. 1, б).

t t + д t t t + д t

а) 6)

Рис. 1. Алгоритмы связывания: а) явный, 6) неявный

Для этого используются, как правило, классические итерационные методы: метод простых итераций с ускорением, метод Гаусса-Зейделя, методы ньютоновского или квазиньютоновского типа и др. [3]. Если придерживаться принятых ограничений, связанных с независимостью от структур данных и методов решения задач в различных физических подобластях, тогда возможно рассмотрение, в основном, только двух методов, которые удовлетворяют этим требованиям, а именно метода Гаусса-Зейделя и квазиньютоновские методы с различными алгоритмами ускорения сходимости.

ПОСТАНОВКА ЗАДАЧИ И УСКОРЕНИЕ СХОДИМОСТИ

Для определенности представления схемы связывания разделенной задачи, будем считать, что на границе-интерфейсе определены нагрузки рг и перемещения иг. Тогда задачу деформирования можно представить через нелинейную функцию S, которая на временном шаге I + Аt, связывает интерфейсные нагрузки рг с перемещениями на интерфейсе 1+миг: ^миг = 8(^рг). Соответственно, газодинамическая задача может быть представлена в виде нелинейной функции F перемещений границы *+Аиг для вычисления нагрузок ^рг = Б^+Аиг).

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

t+Аtuг= 8^+%)] = 8 о Б(^иг). (1)

Для решения нелинейной системы (1) необходимо решить задачи: гидро/газодинамики СБО и динамики деформирования С8О; сопряжения Б81 (обеспечение связи между модулями решения обеих задач, интерполирование давления и деформирование сетки (СМО), реализации семейства итерационных процедур при связывании (рис. 1)).

В предыдущей работе авторов [4] использовался метод Гаусса-Зейделя, который является одним из широко используемых итерационных алгоритмов.

Итак, будем рассматривать решение нелинейной системы для шага по времени I + &t (в дальнейшем временной шаг будет опущен) вида

<+&'иг = 8^+% (*+!) )] = 8 о +% (*+!)). (2)

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

Алгоритм 1. Метод Гаусса-Зейделя с ускорением Эйткена для интерфейсной системы (08)

В обзоре [5] авторы показали, что метод Гаусса-Зейделя может и не сходиться для жесткосвязанных задач Б81, например, при рассмотрении гибких конструкций.

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

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

МОДИФИЦИРОВАННЫЙ МЕТОД БРОЙДЕНА

В настоящее время известно несколько вариантов оригинального метода Бройдена [6]. В дальнейшем рассмотрим комбинацию двух вариантов, сокращающих вычислительные затраты и требуемый объем оперативной памяти. Первый вариант предложен в работе [7]

для задач с ограничением по памяти и использованием перезапуска процесса через задаваемое пользователем максимальное число итераций. Второй - в работе [8], он требует на итерации o(n 2) арифметических операций, а построение аппроксимации якобиана основано на формуле Шермана-Моррисона-Вудбери.

Предлагаемый алгоритм гибридного варианта алгоритма метода Бройдена для решения интерфейсной задачи (2) имеет вид алгоритма 2.

Алгоритм 2. Метод Бройдена c ускорением (BR) 5: В Ь141^1 СЛИ~гьрCaL4M0,4НЬ1 ^ ^^

Рациональным выбором релаксационных параметров W можно добиться существенного ускорения сходимости вычислительного процесса. Способы определения этих параметров основаны на линейной интерполяции, нижней релаксации и использовании асимптотических свойств линейно сходящихся последовательностей [9, 10].

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

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

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

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

Программная реализация решения сопряженной задачи основывается на максимально полном использовании существующего кода, предназначенного для решения отдельных физических задач. Каждая из частей сопряженной задачи решается независимо с использованием гибридных вычислительных узлов. Гидро/газодинамические задачи решаются на основе классов библиотеки ОрепБОАМ. Динамическая задача деформирования решается методом конечных элементов, реализованном в пакете БЕ81иёю [11].

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

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

Рассматривается постановка и алгоритм решения сопряженной задачи взаимодействия сверхзвукового потока и деформируемой панели [1]. Математическая модель взаимодействия сверхзвукового потока газа и деформируемой панели в ударной трубе построена с учетом максимального соответствия условиям проведения эксперимента [12].

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

Деформируемая панель, изготовленная из стали плотностью г = 7600 кг/м и

модулем Юнга Е = 220 ГПа, полагается изотропной и однородной.

Рассмотрены два варианта взаимодействия нестационарной ударной волны с деформируемой панелью в рамках слабосвязанного и жесткосвязанного ^81)

подходов, что дало возможность выполнить непосредственное качественное и количественное сравнение полученных результатов с данными, представленными [13].

На рис. 2 представлены результаты изменения во времени отклонения верхней кромки панели и значений давления в точке «А», в которой регистрировались значения в эксперименте.

Как и в работе [13], смещению кромки против потока соответствуют положительные значения на графиках. Прежде всего, следует отметить хорошее соответствие рассчитанного по представленной модели отклонения панели (рис. 2, а), с экспериментальными данными, так как именно определение деформации конструкции под воздействием потока газа и представляет основной практический интерес в подобного рода задачах. Рассчитанные колебания давления (рис. 2, б), также удовлетворительно согласуются с физической картиной течения, что следует из совпадения положений максимумов, соответствующих прохождению отраженной от панели ударной волны (t = 0,0002 с) и затем волны, отраженной от правой стенки камеры (t = 0,0017 с).

А Па

240 ООО

б)

Рис. 2. Графики зависимостей перемещения края панели (а) и давления (б) от времени: эксперимент [13], расчет wFSI [13], расчет [12], расчет sFSI (алгоритм 1)

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

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

ТЕЧЕНИЕ НЕСЖИМАЕМОЙ ЖИДКОСТИ В КАВЕРНЕ С ГИБКИМ ОСНОВАНИЕМ

Хорошо известная тестовая задача вычислительной гидродинамики была модифицирована для анализа сильносвязанной постановки и тестирования предлагаемого алгоритма Бройдена. Геометрические параметры каверны были следующими: l1 = 1 м, /2 = 0,875 м, l3 = 0,125 м. Постоянная величина горизонтальной скорости набегающего потока была заменена на периодически изменяемую скорость, вычисляемую по соотношению вида uF = 1 — cos (2pt /5). Граничные условия на вертикальных стенках заданы как показано на рис. 3.

Рис. 3. Расчетная область в задаче о течении в каверне с эластичным дном

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

3 2 2

Рр = 1 кг/м и вязкость ьр = 0,01 м /с, а для гибкой мембраны: модуль Юнга Е5 = 250 Н/м , коэффициент Пуассона = 0,0. В расчетах изменялось соотношение плотностей жидкости и деформируемого тела рр / р8 с различными шагами интегрирования по времени, которые

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

Каждый вариант соотношения плотностей считался для промежутка времени ! = 0...50с. Максимально допустимое число внутренних итераций на каждом временном шаге устанавливалось равным 150. Условия сходимости выбирались следующие: '(к+1) / Ыт <е или '+А'Я( к+1) / *+А*~(к+1) —и < е. Результаты расчетов двух алгоритмов

решения нелинейной системы вида (1) для рассматриваемой задачи методом Гаусса-Зейделя и методом Бройдена с различными стратегиями ускорения сходимости сведены в таблицу.

Таблица

Сравнение вычислительных затрат

Рр / Р8 0,005 0,01 0,033 0,005

А! 0,1 0,1 0,05 0,01

ОБ ЦЯ 66 67 67 66 67 66

АТ 6 10 > 150 6 6 5

ВЯ ВЯ 4 6 - 4 4 4

ЦЯ 9 - - 9 9 10

Рассматривались три варианта соотношения плотностей жидкости и эластичной мембраны при трех различных временных шагах интегрирования. В методе Гаусса-Зейделя наряду с ускорением Айткена (см. алгоритм 1) рассматривалось ускорение по методу нижней релаксации ЦЯ. Также метод нижней релаксации использовался в методе Бройдена (см. алгоритм 2) наряду с вариантом без ускорения (ВЯ). В таблице для некоторых расчетных параметров не удалось получить решение (помечено как «-») в связи с вырождением некоторых ячеек вычислительной сетки, что требует дополнительных исследований и пересмотра некоторых алгоритмов ее деформирования.

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

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

Работа выполнена при финансовой поддержке РФФИ (грант 14-01-00055-а, 14-08-00064-а).

СПИСОК ЛИТЕРАТУРЫ

1. Копысов С.П., Тонков Л.Е., Чернова А. А. Двухстороннее связывание при моделировании взаимодействия сверхзвукового потока и деформируемой пластины. Сравнение численных схем и результатов эксперимента // Вычислительная механика сплошных сред. 2013. Т. 6, № 1. С. 78-85.

2. Копысов С.П., Кузьмин И.М., Тонков Л.Е. Алгоритмическое и программное обеспечение решения задач взаимодействия конструкции с жидкостью/газом на гибридных вычислительных системах // Компьютерные исследования и моделирование. 2013. Т. 5, № 2. С. 153-164.

3. Hou G., Wang J., Layton A. Numerical Methods for Fluid-Structure Interaction - A Review // Commun. Comput. Phys. 2012. V. 12, № 2, P. 337-377.

4. Караваев А.С., Копысов С.П., Кузьмин И.М., Недожогин Н.С., Новиков А.К., Рычков В.Н., Сармакеева А.С., Тонков Л.Е., Чернова А.А. Связанные задачи взаимодействия потоков жидкости/газа и деформируемых тел на нестыкующихся сетках // В сб. «Механика и физико-химия гетерогенных сред, наносистем и новых материалов». Ижевск : Изд-во ИМ УрО РАН, 2015. С. 18-41.

5. Degroote J., Bruggeman P., Haelterman R., Vierendeels J. Stability of a coupling technique for partitioned solvers in FSI applications // Computers and Structures. 2008. V. 86. P. 2224-2234.

6. Дэннис Дж., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М. : Мир, 1988. 440 с.

7. Kelley C.T. Iterative methods for linear and nonlinear equations. Philadelphia : SIAM, 1995. 166 p.

8. Kelley C.T. Solving nonlinear equations with Newton's method. Philadelphia : SIAM, 2003. 105 p.

9. Трауб Д.Ф., Вожьняковский Х. Общая теория оптимальных алгоритмов. М. : Мир, 1983. 382 с.

10. Вержбицкий В.М. Основы численных методов. М. : Высшая школа, 2002. 840 с.

11. Копысов С.П., Кузьмин И.М., Недожогин Н.С., Новиков А.К., Рычков В.Н., Сагдеева Ю.А., Тонков Л.Е. Параллельная реализация конечно-элементных алгоритмов на графических ускорителях в программном комплексе FEStudio // Компьютерные исследования и моделирование. 2014. Т. 6, № 1. С. 79-97.

12. Копысов С.П., Кузьмин И.М., Тонков Л.Е. Моделирование взаимодействия сверхзвукового потока и деформируемой панели в ударной трубе // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2012. № 2. С. 156-165.

13. Giordano J., Jourdan G., Burtschell Y., Medale M., Zeitoun D.E., Houas L. Shock wave impacts on deforming panel, an application of fluid-structure interaction // Shock Waves. 2005. V. 14. P. 103-110.

QUASI-NEWTON METHOD FOR IMPLICIT COUPLING IN FSI PROBLEM

Kopysov S.P., Kuzmin I.M., Nedozhogin N.S., Novikov A.K., Rychkov V.N., Tonkov L.E. Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia

SUMMARY. We consider split approach for solving strong coupling fluid-structure interaction problem on the non-matching meshes. We offer algorithm for solve nonlinear equations system on the interaction border. We present modeling results within strong coupling split approach of the some test problem.

KEYWORDS: coupling problems, nonlinear system of equation, Broiden method, convergence acceleration, fluid-structure interaction problem.

Копысов Сергей Петрович, доктор физико-математических наук, заведующий лабораторией ИМ УрО РАН, e-mail: s. kopysov@gmail. com

Кузьмин Игорь Михайлович, младший научный сотрудник ИМ УрО РАН, e-mail: imkuzmin@gmail com

Недожогин Никита Сергеевич, научный сотрудник ИМ УрО РАН, e-mail: nedozhogin@inbox.ru

Новиков Александр Константинович, кандидат физико-математических наук, старший научный сотрудник ИМ УрО РАН, e-mail: sc_work@mail. ru

Рычков Владимир Николаевич, кандидат физико-математических наук, старший научный сотрудник ИМ УрО РАН, e-mail: bob.r@mail.ru

Тонков Леонид Евгеньевич кандидат физико-математических наук, старший научный сотрудник ИМ УрО РАН, e-mail: tnk@@udman. ru