CC BY
38
УДК 534.141,533.6
М. Л. Дмитриев, Г. А. Щеглов
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ УПРУГОЙ СИСТЕМЫ ПРИ ЕЕ ОТРЫВНОМ ОБТЕКАНИИ ИДЕАЛЬНОЙ НЕСЖИМАЕМОЙ СРЕДОЙ
Предложена модель динамики осциллятора, состоящая из подпружиненного тела, обтекаемого плоскопараллельным потоком идеальной жидкости. Рассмотрено нестационарное плоскопараллельное отрывное обтекание тела. В качестве метода моделирования использован метод дискретных вихрей. Приведены некоторые результаты моделирования в сравнении с результатами экспериментов, выполненных другими авторами.
При решении задач гидроупругости возникает необходимость численного моделирования явлений, обусловленных отрывом потока с поверхности обтекаемых тел. При определении мест отрыва потока традиционными являются подходы, связанные с учетом вязкости среды, а именно: решение уравнений Навье-Стокса во всей рассматриваемой области течения или/и решение уравнений пограничного слоя вблизи поверхности обтекаемого тела [1, 2]. При больших числах Рейнольдса (106 и более) используются модели турбулентности, имеющие, как правило, полуэмпирический и/или неуниверсальный характер. Кроме того, для проведения исследований необходимы значительные затраты времени счета и существенные вычислительные возможности ЭВМ.
Для эффективного решения инженерных задач на практике интерес также представляют подходы к моделированию отрывных течений с использованием более простых моделей среды, например модели обтекания тел в рамках теории идеальной жидкости. Особый интерес к таким методам возникает в случае решения задач гидроупругости, когда одновременно с определением аэродинамических нагрузок на конструкцию необходимы динамические расчеты, также требующие значительных вычислительных ресурсов. В настоящей работе в целях определения гидродинамических нагрузок использован нашедший широкое применение метод дискретных вихрей [3].
Рассмотрим в общем случае неустановившийся режим движения упруго закрепленного твердого тела в идеальной несжимаемой среде. В качестве модельной выбрана задача о колебаниях системы с одной степенью свободы, состоящей из кругового цилиндра массой М, закрепленного на пружине жесткости С (рис. 1). В начальный момент времени колебательная система покоится в положении равновесия в
Рис. 1. Расчетная схема задачи
набегающем потоке, имеющем постоянную скорость. Отрывное обтекание моделируется с помощью вихревых пелен, сходящих с поверхности тела в каждый момент времени в точках отрыва (точки А и В на рис. 1). Течение вне поверхности тела и его вихревого следа считается потенциальным. Граничные условия: непротекание на поверхности тела, совмест-
ность течения на поверхностях вихревых пелен и затухание возмущений на бесконечном удалении от тела и его следа. Поле давлений в течении вычисляется на основе интеграла Коши-Лагранжа [2].
Упругая подсистема описывается уравнением малых колебаний гармонического осциллятора
где у — вертикальное перемещение центра масс цилиндра; и2 = С/М; Р (Ь,у,у) — проекция нестационарной гидродинамической силы Р на ось ОУ. Уравнение (1) при моделировании решается методом Рунге-Кутты четвертого порядка.
При численном определении гидродинамической силы методом дискретных вихрей поверхности тела и пелен моделируется система точечных вихрей. На поверхности тела в контрольных точках, находящихся между точечными вихрями, выполняется условие непротекания, приводящее в каждый момент времени к системе линейных алгебраических уравнений относительно неизвестных интенсивностей вихрей 7, моделирующих тело:
где [А] — матрица с элементами а^ = щ • п, щ — скорость, индуцированная вихрем единичной интенстивности, находящимся в точке с номером ], в контрольной токе с номером г; V — скорость частиц жидкости на поверхности тела в неподвижной системе координат ОХУ (см. рис. 1); г* — радиус-вектор точки тела; п — внешняя нормаль к поверхности тела.
Интенсивности точечных вихрей, пополняющих в каждый момент времени систему вихрей пелен, определяются из условия схода в поток всей циркуляции вихревого слоя в точках отрыва. Эволюция вихревых пелен описывается системой дифференциальных уравнений движения
y + и2y = F (t,y,y),
(1)
точечных вихреи по траекториям жидких частиц, которая численно решается методом ЭИлера;
dr ->
dt = V ^ (2)
Гидродинамическая сила, действующая на осциллятор, определяется в результате интегрирования сил давления по поверхности тела в каждый момент времени по методике, используемой авторами работы [2] в случае плоскопараллельного течения.
Местоположение точек отрыва в каждый момент времени находится на основе предложенного В.А. Апариновым и М.Л. Дмитриевым [4, 5] критерия о минимальности затрат энергии на создание нового элементарного участка вихревой пелены, сходящей в поток в возможной точке отрыва.
В рассматриваемый момент времени вихревые пелены, с помощью которых моделируется отрывное обтекание, сходят с поверхности в точках с криволинейными координатами sa, sb (см. рис. 1), где достигаются минимальные значения:
K (s) = V(s)2(v(s)„n(s)).
—t —t dr
Здесь V = V--— относительная скорость частиц жидкости на
dt ^
поверхности тела; v = ~г — Voo, где V^ — скорость набегающего по-
dt
тока в неподвижной системе координат ОXY (рис. 1). Критерий K (s) получен вариацией уравнения энергии и характеризует работу, затраченную в единицу времени на создание возможного элементарного участка вихревой пелены в точке с координатой s [4, 5]. В данной задаче использована расчетная схема с двумя вихревыми пеленами. Следовательно, в каждый момент времени имеются две точки отрыва sa,sb , соответствующие двум первым глобальным минимумам на контуре цилиндра.
Исследование системы (1) включало в себя два этапа изучения динамических характеристик при различных видах колебаний тела в потоке.
На первом этапе считали, что C = 0 и M = 0.В целях тестирования методики было проведено моделирование режимов нестационарного отрывного обтекания цилиндра, совершающего перемещения поперек набегающего потока по гармоническому закону с амплитудой, равной одной четверти диаметра цилиндра:
y(t) = 0,25 d cos ш t.
Полученные в ходе численных экспериментов результаты сравнивались с известными экспериментальными данными [6], показываю-
щими, что имеется зависимость между частотой ш колебаний цилиндра и частотой f изменения давления в заданной точке поверхности цилиндра. Закон движения профиля в потоке невязкой несжимаемой среды соответствовал описанным условиям натурного эксперимента [6] с учетом перехода к безразмерным параметрам (V» = 1, р^ = 1, р =1, d = 1). На рис.2 показаны отдельные фазы колебаний профиля в потоке и соответствующие им фазы формирования вихревой пелены в различные моменты безразмерного времени. Также на рисунке показаны эпюры критерия К (з), положения точек отрыва и главный вектор гидродинамических сил 1.
Численные эксперименты проводились для частот ш, соответствующих диапазону чисел Струхаля 0 < ш < 0,3. В ходе вычислений строились зависимости от времени коэффициента давления в точках с координатами ф = ±30°; ±40, 5°, расположенных на поверхности цилиндра. Пример зависимости коэффициента давления от времени для ф = +40,5° показан на рис. 3 (время указано безразмерное).
Рис. 2. Эволюция гидроупругой системы
CPd
10
Рис. 3. Осциллограмма пульсаций коэффициента давления на поверхности цилиндра в точке ф = 40,5° при = 0,1
Рис. 4. Зависимость частоты пульсаций давления 8И$ на поверхности цилиндра от частоты колебаний цилиндра :
сплошные линии — эксперимент [6]; точки — результаты численного моделирования
Sh„ =
В результате анализа спектра пульсаций давления для каждого эксперимента выявлено число БИ^, соответствующее числу БИ^. Сопоставление полученных результатов с данными натурных экспериментов (сплошная кривая) приведено на рис.4.
В эксперименте четко выделяются три качественно различных режима пульсаций давлений в зависимости от числа Струхаля
У^'
где ( — диаметр цилиндра; Vэо — скорость набегающего потока.
Первый режим реализуется при малой частоте колебаний цилиндра (0 < БИ^ < 0,04). При этом частота пульсаций давления, выражаемая числом Струхаля БИ^, не зависит от частоты колебаний цилиндра. Второй режим (0,04 < БИ^ < 0,1) характеризуется снижением частоты пульсаций давления при увеличении частоты колебаний цилиндра. Третий режим (БИ^ > 0,1) характеризуется синхронизацией частоты колебаний давления с частотой колебаний цилиндра. Полученная в
эксперименте зависимость числа Струхаля БИ/ = fd/V(Xl для пульсаций давления от числа БЪ^ приведена на рис. 4 сплошной кривой.
Сравнение данных, полученных численно, с экспериментальной зависимостью показывает, что все три качественно различных режима отрывного течения верно моделируются с использованием энергетического критерия нахождения точки отрыва. Полученные результаты свидетельствуют о возможности применения рассматриваемой математической модели для исследования характеристик неустановившегося движения гидроупругой динамической системы.
На втором этапе предлагаемая в настоящей работе методика была применена непосредственно к исследованию колебаний осциллятора (см. рис. 1), возникающих под действием силового возбуждения потоком среды, имеющим единичные безразмерные параметры ^^ = 1, р^ = 1, р = 1). При этом проекция вектора силы 1 (см. рис.2) на ось ОУ была подставлена в правую часть уравнения (1). В качестве примера далее приведены результаты расчета переходных режимов для круглого профиля диаметром d = 1,0; массой М = 3,0 (величины безразмерные).
На рис. 5, а показан переходной режим движения центра масс упруго закрепленного профиля в случае, когда пружина имеет безразмерную жесткость С = 2,0. Собственная частота колебаний такой системы БЪ = 0,13 ниже частоты схода вихрей (БЪ = 0,20) и амплитуда вынужденных колебаний невелика.
На рис.5,б показан переходной режим в случае, когда пружина имеет безразмерную жесткость С = 4,737. Собственная частота колебаний системы (БЪ = 0,2) практически совпадает с частотой схода вихрей и амплитуда вынужденных колебаний становится большой. Однако неограниченного роста амплитуды не происходит из-за нелинейного взаимодействия профиля с потоком.
В системе можно также наблюдать режимы биений, как показано, например, на рис.5,в для случая С = 20 (собственная частота осциллятора соответствует БЪ = 0,41).
С увеличением жесткости пружины амплитуда вынужденных колебаний уменьшается. На рис.5,г показан график движения центра профиля для случая С = 118,43 (собственная частота осциллятора соответствует БЪ = 1,00).
На примере проведенных расчетов видно, что методика позволяет моделировать основные эффекты взаимодействия упругой конструкции с потоком.
Основным выводом настоящей работы следует считать успешное тестирование методики определения гидродинамических нагрузок на колебательную систему с использованием метода дискретных вихрей и энергетического критерия определения мест отрыва потока. Данная
Рис. 5. График движения осциллятора с собственной частотой, соответствующей 8И = 0,13 (а); 0,20 (б); 0,41 (в) и 1,00 (г)
методика позволяет моделировать переходные режимы с учетом взаимодействия упругой конструкции с потоком идеальной несжимаемой среды.
Существенным является также тот факт, что машинное время, необходимое для проведения 2050 шагов интегрирования системы уравнений (1) и (2) при расчете каждого варианта на персональном компьютере с процессором Pentium IV (2 ГГц) не превышало 2 ч. Такие временные затраты на порядок меньше, чем при использовании сеточных методов, что характеризует рассмотренный метод как весьма эффективный.
Предложенный в настоящей работе подход предполагается использовать в дальнейшем при исследовании динамики более сложных двумерных и трехмерных аэроупругих систем.
СПИСОК ЛИТЕРАТУРЫ
1. Р о у ч П. Вычислительная гидродинамика. - М.: Мир, 1980. - 616 с.
2. Трехмерное отрывное обтекание тел произвольной формы / Под ред. С.М. Белоцерковского. - М.:ЦАГИ, 2000. - 265 с.
3. Нелинейная теория крыла и ее приложения // Т.О.Аубакиров, С.М. Белоцерковский, А.И. Желанников, М.И. Ништ. - Алматы: Гылым, 1997. - 448 с.
4. Д м и т р и е в М. Л. Математическое моделирование плоскопараллельного отрывного обтекания тел потоком идеальной жидкости // НММ "Аэродинамика летательных аппаратов". - М.: ВВИА им. Н.Е. Жуковского, 1996.
5. С п у т н ы е следы и их воздействие на летательные аппараты. Моделирование на ЭВМ / Т.О. Аубакиров, А.И. Желанников, П.Е. Иванов, М.И. Ништ. -Алматы: Гылым, 1999. - 280 с.
6. Федяевский К. К., Блюмина Л. X. Гидродинамика отрывного обтекания тел. - М.: Машиностроение, 1977.
Статья поступила в редакцию 31.01.2007
Михаил Леонардович Дмитриев родился в 1968 г., окончил Московский физико-технический институт в 1991 г. Канд. техн. наук, главный специалист ООО "УНИИ АРКС". Автор 20 научных работ в области численного моделирования в аэрогидродинамике.
M.L. Dmitriev (b. 1968) graduated from the Moscow Physical Technical Institute in 1991. Ph. D. (Eng.), chief specialist of the limited stock company "OOO "UNII ARKS". Author of 20 publications in the field of numerical simulation in aero-gasdynamics.
Георгий Александрович Щеглов родился в 1972 г., окончил МГТУ им. Н.Э. Баумана в 1996 г. Канд. физ.-мат. наук, доцент кафедры "Аэрокосмические системы" МГТУ им. Н.Э. Баумана. Автор 28 научных работ в области динамики и прочности космических летательных аппаратов и прикладной механики.
G.A. Shcheglov (b. 1972) graduated from the Bauman Moscow State Technical University in 1996. Ph. D. (Phys.-Math.), assoc. professor of "Aerospace Systems" department of the Bauman Moscow State Technical University. Author of 28 publications in the field of dynamics and strength of space flying vehicles and applied mechanics.
