Задачи вычислительной гидродинамики (CFD) относят к наиболее сложным задачам вычислительной математики в силу нелинейности определяющих соотношений. Если в инженерных приложениях возникают сопряженные задачи в гидроупругой постановке, когда обтекаемая конструкция взаимодействует с потоком, двигаясь или деформируясь под действием аэрогидродинамических нагрузок (Fluid-Structure Interaction, FSI), задача усложняется еще сильнее, поскольку форма области течения становится переменной, а ее границы – подвижными, причем их скорость в общем случае заранее неизвестна. В настоящей работе рассмотрены алгоритмы решения таких задач с использованием вихревых методов вычислительной гидродинамики [1], которые могут быть чрезвычайно эффективны, если влиянием сжимаемости среды можно прененбречь, а подвижный элемент конструкции имеет значительное удлинение, так что его обтекание можно считать плоским.
Вихревые методы относятся к бессеточным лагранжевым методам моделирования течений, в которых первичной расчетной величиной является завихренность . Моделируя эволюцию (перенос) завихренности в области течения и ее генерацию на обтекаемой поверхности, можно описать течение вязкой несжимаемой жидкости вокруг профиля.
Вихревые методы в силу своей низкой вычислительной сложности хорошо применимы для моделирования существенно нестационарных режимов обтекания профилей или их систем (в том числе подвижных или деформируемых). При этом моделирование обтекания с применением вихревых методов позволяет решать задачу в наиболее общей постановке, отказавшись от гипотезы (квази)стационарности нагрузок, т. е. их вычисления по заранее найденным значениям стационарных гидродинамических коэффициентов, и рассматривая непотенциальный вихревой режим обтекания.
Благодаря лагранжевой природе вихревых методов и отсутствию необходимости построения расчетной сетки в области течения вихревые методы эффективны в задачах, где форма области течения изменяется во времени; при этом обтекаемая поверхность может совершать произвольные движения или деформироваться.
Сложность представляет организация схемы связывания упругой и гидродинамической подсистем; фундаментальные результаты, полученные в [2], устанавливают связь между величинами гидродинамических нагрузок и интенсивностью генерируемой на обтекаемой поверхности завихренности, которая оказывается линейной. Это позволяет построить полунеявную разностную схему интегрирования по времени, итерации в которой сходятся достаточно быстро, что открывает путь к устойчивому и эффективному решению задач для тел произвольной массы, вплоть исчезающе малой.
На сегодняшний день разработан программный комплекс VM2D [3] для решения плоских задач моделирования несжимаемых течений и расчета обтекания профилей. В его основу положен метод вязких вихревых доменов [4], а также некоторые оригинальные разработки авторов, связанные, главным образом, со схемами численного решения граничных интегральных уравнений относительно интенсивности генерируемого на профилях вихревого слоя [5], рис. 1.
На текущий момент комплекс VM2D позволяет решать широкий класс задач расчета обтекания подвижных и неподвижных профилей и систем профилей и определения действующих на них нестационарных гидродинамических нагрузок, однако решение сопряженных задач гидроупругости в случаях, когда инерционные свойства (масса, момент инерции) профиля соизмеримы с таковыми для вытесняемого объема среды (такие профили естественно называть «легкими»), приводит к значительной погрешности или даже вычислительной неустойчивости. Это является следствием применения наиболее простой схемы расщепления, в соответствии с которой каждый шаг расчета по времени делится на два подшага: на первом рассчитывается обтекание профиля, движущегося с известной скоростью (известной с предыдущего шага), на втором — движение профиля под действием известных гидродинамических нагрузок (найденных на первом подшаге).
В рамках настоящей работы в пакет VM2D внедрен итерационный алгоритм моделирования гидроупругих колебаний профилей. Данный алгоритм основан на результатах [2], и требует для реализации полунеявной схемы расчета компонент тензора присоединенных масс профиля, которые могут быть легко вычислены с использованием вышеупомянутых схем решения граничных интегральных уравнений [5].
Работа реализованного итерационного алгоритма продемонстрирована на модельной задаче о расчете гидроупругих колебаний, в которой моделируется обтекание потоком вязкой несжимаемой жидкости кругового профиля, закрепленного посредством вязкоупругой связи и способного перемещаться поперек потока. Частота собственных колебаний механической системы выбрана близкой к частоте схода вихрей Кармана, что соответствует режиму ветрового резонанса.
Зависимость отклонения «тяжелого» профиля от времени показана на рис. 2; результаты моделированияблизки для ранее реализованного в VM2D и обсуждаемого здесь алгоритмов.
Для профиля, плотность которого близка к плотности вытесняемой им среды (рис. 3), амплитуда колебаний при расчете обоими алгоритмами оказывается близкой, однако частота колебаний отличается весьма значительно, причем новый алгоритм обеспечивает более низкую частоту колебаний, что качественно лучше соответствует поведению такой механической системы.
Для моделирования колебаний «легкого» профиля, масса которого меньше массы вытесняемой им среды (а точее — его присоединенной массы), ранее реализованный в VM2D алгогитм был непригоден, посколько приводил к вычислительной неустойчивости, тогда как новый алгоритм позволяет корректно решать задачу, требуя, в зависимости от конкретной величины массы профиля, от единиц до нескольких десятков итераций до сходимости на каждом шаге расчета по времени. Отметим, что в случае расчета обтекания «тяжелых» профилей требуется, как правило, не более двух таких итераций, а для экономии времени расчета можно ограничиться и безытерационной схемой.
Реализованный в программном комплексе VM2D алгоритм [2] обеспечивает возможность устойчивого моделирования в сопряженных задачах гидроупругости при произвольной массе/моменте инерции обтекаемого профиля и возможность получения качественно правильного результата, в частности, зависимости частоты колебаний от относительной массы профиля. Полунеявный характер схемы расщепления позволяет проводить расчет с достаточно крупным шагом по времени и требует выполнения небольшого количества итераций до сходимости соответствующего процесса.