CC BY
80
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
Сер. 10. 2009. Вып. 3
УДК 629.12.035
Д. В. Никущенко, А. Л. Чистов
АЛГОРИТМ РАСЧЕТА ГИДРОДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК СИСТЕМ КРЫЛЬЕВ ОТНОСИТЕЛЬНО БОЛЬШОЙ ТОЛЩИНЫ
При решении задач динамики морских подводных объектов часто возникает необходимость производить расчеты гидродинамических характеристик систем крыльев (к числу последних относятся прежде всего органы управления, стабилизации и ограждения выдвижных устройств). Главные особенности таких крыльев - относительно большая толщина профиля и малое удлинение.
Одним из популярных методов расчета крыльевых систем является метод дискретных вихрей (МДВ) - разновидность метода гидродинамических особенностей. Большой вклад в его разработку внесли В. А. Апаринов, С. М. Белоцерковский, А. А. Зайцев, А. В. Дворак, И. К. Лифанов, М. И. Ништ, Н. Н. Поляхов, В. К. Трешков, F. S. Archibald, H. Ashley, P. T. Fink, W. Kraus, A. Leonard, R. I. Lewis, P. Koumoutsakos, P. G. Saffman, T. Sarpkaya и др. Этот метод показал очень хорошие результаты для крыльев самолетов, однако расчет крыльев, применяемых на подводных объектах, имеет ряд особенностей [1]. В частности, малое удлинение таких крыльев (как правило, 1-1.5) приводит к тому, что между потоками, сходящими с торцевых кромок крыла, возникает сильное гидродинамическое взаимодействие, оказывающее сильное влияние на его гидродинамические характеристики. В то же время вследствие большой относительной толщины (18-21% хорды профиля В) происходит изменение эпюры распределения давлений по поверхности крыла в сравнении с авиационными крыльями. Большинство же существующих алгоритмов расчета крыльевых систем предназначено для моделирования обтекания крыльев конечного удлинения.
В настоящей работе представляется расчетный комплекс, предназначенный для расчета гидродинамических характеристик крыльевых систем морских подводных объектов. Численная схема, примененная для создания данного комплекса, позволяет учитывать указанные особенности обтекания крыльев, используемых на подводных объектах.
При решении поставленной задачи жидкость считается идеальной и невесомой, а движение ее вне вихревых поверхностей - потенциальным: v = УФ, где Ф(г, t) - потенциал вызванных скоростей в данной точке пространства.
Поверхность каждого моделируемого объекта заменяется непрерывным распределением вихревых особенностей (присоединенных вихрей), тогда в поток с нее должны
Никущенко Дмитрий Владимирович — доцент кафедры теории корабля факультета кораблестроения и океанотехники Санкт-Петербургского государственного морского технического университета. Количество опубликованных работ: 46. Научные направления: теория турбулентности, динамика судов. E-mail: [email protected].
Чистов Алексей Леонидович — аспирант кафедры вычислительных методов механики деформируемого тела факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Научный руководитель: проф. В. А. Павловский. Количество опубликованных работ: 3. Научные направления: гидродинамика, теория турбулентности. E-mail: [email protected].
© Д. В. Никущенко, А. Л. Чистов, 2009
сходить свободные вихри, обусловленные изменением интенсивности присоединенных (т. е. местной скорости потока).
В соответствии с уравнением неразрывности V- v = У-УФ = 0, вне вихревых поверхностей (т. е. поверхности тела 5 и вихревой пелены Бр (рис. 1)) потенциал скоростей Ф является гармонической функцией (т. е. удовлетворяет уравнению Лапласа).
Рис. 1. Вихревая система прямоугольного крыла
Таким образом, задача об обтекании тела потоком невязкой жидкости сводится к отысканию потенциала вызванных скоростей Ф(г, і), удовлетворяющего уравнению Лапласа
ДФ = 0
при следующих граничных условиях:
1) на поверхности объекта Я выполняется условие непротекания
д Ф(і)
д п(і)
= ^ (і) • п(і),
Б
где (Ь) - скорость движения точек поверхности 5 (или скорость набегающего потока
при обращенном движении);
2) на бесконечности выполняется условие убывания возмущений
|VФ (гм^)| ^ 0 при гм ^ го;
3) на пелене свободных вихрей Бр осуществляются условия совместности течения
тйг = р+ <м»> = »■<"»>• м»е^
4) на линии схода вихревой пелены исполняется условие Чаплыгина-Жуковского о конечности скоростей (т. е. градиент Ф ограничен в окрестностях точек Мо Е Б в любой момент времени Ь).
В результате выполнения условий непротекания и убывания возмущений задача сводится к решению интегродифференциального уравнения [2]
7-7Г— Е / 7Г~ (~~= _у°° ' пмА*)> (!)
4п дпм0^ ] дпм \тмм0)
г=1 Я
где N - общее количество вихревых поверхностей; - интенсивности вихревых слоев; М0 - точка, принадлежащая вихревой поверхности; М - точка вне вихревых поверхностей.
Детальное рассмотрение вопросов существования и единственности решения приведенного уравнения содержится в монографии И. К. Лифанова [2].
Окончательно задача сводится к решению уравнения (1) относительно неизвестных интенсивностей присоединенных вихревых особенностей ^г, по которым возможно определять интенсивности свободных вихрей, если геометрия вихревых пелен известна. Отсутствие нагрузки на пелене свободных вихрей (краевое условие 3) приводит к следующему: точки вихревой пелены должны перемещаться по направлению местной скорости. То есть для точки Мр(Ь,Ьо), сошедшей с линии Ьр в момент Ьо, в момент £ должно выполняться следующее равенство:
= У ^Мр ^ ^ + ^Мр ^ ^ + У°° ^Мр ^ ^ ’ где V - вызванная скорость в точке Мр от системы присоединенных вихрей; vp - вызванная скорость в точке Мр от системы свободных вихрей на временном слое Ь.
Поставленная задача может решаться в связанной системе координат для случая как обращенного движения (т. е. когда неподвижное тело обтекается потоком жидкости), так и абсолютного движения (т. е. когда тело движется с поступательной скоростью v0(t) и угловой ш(Ь) в изначально неподвижной жидкости).
Для решения уравнения (1) непрерывное распределение вихрей заменяется дискретным, для чего поверхности крыльев системы продольными и поперечными линиями разбиваются на элементарные площадки - панели (это можно видеть на рис. 1). Воздействие каждой панели на поток заменяется воздействием вихревой рамки постоянной интенсивности. Тогда уравнение (1) заменяется системой алгебраических уравнений относительно интенсивностей вихревых рамок:
N т №р к
+^^JГjaqj = —4п"МпІ — ЕЕЕ 5иакрГ+1, з = I, ч = 1, 2,...,М, (2)
= 1 р=1 І=1 в=1
в которой wnj = vnj/ух - безразмерная нормальная составляющая скорости невозмущенного потока в точке Mj, ГІ = Г+/(|уто|В) - безразмерная циркуляция г-й вихревой рамки, ХN - регуляризирующая переменная; I - номер рамки с нулевой циркуляцией, 5р)і = Гт-1 — Гк-1 - циркуляция рамки, сходящей с линии Ьр (равная разности циркуляций прилегающих рамок), а^ - координатная функция, определяемая взаимным расположением точек Mj и МІ на поверхности тела 5:
4 . .
___________{гк+і-гк)-{гк—г:і)__________ _ І (гк+і-гк)-(гк-гі) _ (гк+і—гк)-{гк—г;і) 1
(гк-^)2(гк+1-Гк)2-[(гк+1-Гк)(гк-г^]2 \ \гк+1-Гз\ \гк-гз\ / '
а
Вихревую рамку с нулевой циркуляцией необходимо вводить в силу того, что в противном случае СЛАУ становится вырожденной [3]. Получившаяся в итоге система (2) является переопределенной, поэтому в работе [4] было предложено вводить регуляри-зирующую переменную XN.
После вычисления значений циркуляций вихревых рамок можно рассчитать скорости потока во всей области течения. Затем, используя интеграл Коши-Лагранжа или теорему Н. Е. Жуковского «в малом»[5], определить гидродинамические силы во всех элементарных площадках поверхностей тел системы, после чего найти силы и моменты, действующие со стороны жидкости как на всю систему целиком, так и на каждое тело в составе системы.
По изложенному алгоритму Д. В. Никущенко совместно с аспиранткой Е. А. Рогожиной был создан программный расчетный комплекс WingSim (Wings’ Simulation) [1, 6-8]. Он позволяет выполнять моделирование как относительного, так и абсолютного движения системы крыльев, причем движение может быть и поступательным, и вращательным. Возможно использование различных комбинаций вихревых пелен, сходящих с ее элементов - только с задних кромок («линеаризованный» алгоритм метода дискретных вихрей), а также с торцевых и передней (нелинейный алгоритм). Комплекс основан на концепции объектно-ориентированного программирования и реализован на языке программирования Object Pascal в среде программирования Delphi (рис. 2).
111ШШ!111ШШ1!111111111111111111111111111111111111111111111Д И Шщ
| * Вычисление циркуляции сходящих рамок &
| * Расчет поля скоростей следа
| * Определение нового положения вихревого след»
1НАГМ
| Расчет Силы и момента элемента системы:1
| Вычисление результирующих
1* Определение правых частей СЛАУ
1* Решение СЛАУ разложением матрицы по
гашгулярнъхм числам (SV.Il>)
1* Вычисление циркуляции сходящих рамок
I * Распет поля скоростей следа
I* Определение нового положения вихревого следа
| Создан файл геометрии следаС:\Тетр^$.са$1.$Ь(
Рисование тел системы
Расчет завершен
Рис. 2. Главное окно расчетного комплекса WingSim
В соответствии с такой концепцией для построения рассчитываемых тел вводится в рассмотрение базовый класс, которым в данном случае является класс СРапе1, представляющий собой собственно одну вихревую рамку и содержащий, с одной стороны,
ряд полей, задающих положение ее узловых и контрольной точек, векторы единичной нормали и вызванной скорости, гидродинамические силу и момент, действующие на рамку, и, с другой - методы, позволяющие осуществлять операции с указанными полями данных (прежде всего вычисление вектора нормали и вызванной скорости в заданной точке). Каждый моделируемый объект представляет собой класс, который является потомком класса СРапе1, в частности, с точки зрения программы отдельное крыло - экземпляр класса CWing, содержащий большой набор полей и методов, дающих возможность создавать и производить расчеты обтекания крыльевого элемента.
Рассчитываемая система представляет собой экземпляр класса Т8уБ1еш, включает в себя методы, осуществляющие расчеты геометрических параметров системы, создание и решение матрицы вызванных скоростей, вычисление распределенных и интегральных параметров системы в целом и каждого ее элемента в отдельности, и может содержать любое количество экземпляров класса CWing, т. е. крыльев, или других классов, реализованных в программе (рис. 3). Задание геометрических параметров элемента рассчитываемой системы выполняется в локальной системе координат (рис. 4). Кроме того, необходимо задать положение начала локальной системы элемента в глобальной системе (в качестве которой принята обычная система, используемая для расчетов динамики подводных объектов: ось х - вперед, ось у - вверх, ось г - на правый борт). Программа также позволяет задавать пары элементов, подобных приведенным на рис. 4.
Рис. 3. Окно задания элементов рассчитываемой системы
После создания моделируемой системы возможно ее сохранение в файле специального вида вместе с иными параметрами расчета, например угла атаки, скорости, расчетного времени. Комплекс WingSim содержит различные способы решения системы
Рис. 4. Окно задания геометрических параметров составного элемента
алгебраических уравнений метода дискретных вихревых особенностей (МДВО) (рис. 5, А). Это позволяет выполнять расчеты сложных систем, используя наиболее подходящий подход, от метода Гаусса до итерационных [9].
На рис. 5, Б показано окно выбора способа расчета интегральных характеристик системы, выбора точки привязки моментов гидродинамических сил. После выполнения расчета можно просмотреть вихревую систему, получившуюся в ходе моделирования, как для всей системы в целом, так и для ее элементов по отдельности (рис. 6, А).
Таким образом, процедуру решения задачи обтекания системы объектов можно представить в виде следующих шагов:
1) создание расчетной системы в соответствующем окне (см. рис. 2);
2) после нажатия кнопки «Принять» происходит инициализация экземпляров созданных объектов (выделение памяти для полей данных), на этом этапе также рассчитываются параметры вихревых систем - положения узловых и контрольных
Рис. 5. Окно выбора способа решения системы уравнений (А) и способа интегральных характеристик (Б)
точек, внешних нормалей, положения линий схода свободных вихрей;
3) задание параметров расчета - угла атаки, расчетного времени, способов вычисления сил и решения СЛАУ;
4) формирование матрицы координатных функций а^;
Рис. 6. Коэффициенты нормальной силы (А) и продольного момента (Б) для крыла удлинением 1.5 по результатам расчета (1), по данным справочника [10] (2) и работы [11] (3 - линейная вихревая модель, 4 - нелинейная модель )
5) решение системы (2) на нулевом шаге (бесциркуляционной задачи) относительно неизвестных циркуляций вихревых отрезков;
6) определение координат вихревых рамок, сошедших с линий схода Lp, и их интенсивностей 5 + pik;
7) определение суммарных вызванных скоростей в контрольных точках (с учетом влияния вихревой пелены);
8) циклическая процедура решения системы (2) до достижения заданного количества расчетных шагов;
9) вычисление гидродинамических характеристик системы с помощью интеграла Коши-Лагранжа или теоремы Н. Е. Жуковского «в малом»;
10) вывод вида вихревой системы и значений гидродинамических характеристик по командам пользователя.
В целом можно сказать, что созданный одним из авторов расчетный комплекс WingSim представляет собой эффективное средство моделирования обтекания систем разнородных элементов, таких как крылья или замкнутые тела произвольной формы (для этого в нем реализован алгоритм импорта расчетных сеток сеточного генератора «Gambit ®»). Комплекс легко может модернизироваться путем введения новых классов, реализующих операции с дополнительными элементами, как новыми, так и наследуемыми от существующих (например, класса CWing). В текущей редакции комплекс дает возможность также проводить расчет силы волнового сопротивления системы, движущейся вблизи границы раздела сред.
Данный алгоритм позволяет производить расчеты крыльев при различных комбинациях пелен свободных вихрей - пелена может сходить и с задней (этот вариант в алгоритме осуществляется всегда), и с передней кромок, с одной стороны, и с каждой из боковых кромок или с передней - с другой.
Процедура моделирования с помощью комплекса WingSim иллюстрируется на примере обтекания крыла с относительным удлинением А = 1.5, имеющего профиль NACA-0006. Поверхность крыла разбивалась на 320 панелей. Результаты моделирования приведены на рис. 6, 7.
Рис. 7. Установление интегральных характеристик Су (1) и шг (2) в зависимости от номера итерации (расчетного времени)
Следует отметить хорошее согласование коэффициентов нормальной силы (рис. 6, А). Некоторое расхождение коэффициентов момента (рис. 6, Б) можно объяснить тем, что использованный алгоритм не позволяет производить расчеты крыльев предельно малой толщины (относительной толщины f = t/b ^ 0, результаты расчетов которого приведены в справочнике [10]), так как вихревые элементы располагаются по поверхности объекта и при уменьшении толщины крыла вихревые отрезки оказываются слишком близко друг к другу. Как показывают полученные результаты, наименьшее значение относительной толщины при применении данного подхода составляет t/В = 0.06.
На рис. 7 можно видеть установление интегральных характеристик в зависимости от расчетного времени (номера итерации). Для получения решения при моделировании обтекания изолированного крыла необходимо выполнять не менее 60-65 итераций, т. е. производить расчет на интервале расчетного времени не менее 4.5 с для крыла хордой 1 м, движущегося со скоростью 1 м/с.
Время расчета одной точки при разбивке крыла на 320 панелей на 80 расчетных шагах на компьютере AMD Athlon X2 Barton 3800+, 1 Gb ОЗУ составляло в среднем 9 мин в линейном приближении и 38 мин в нелинейном.
Таким образом, можно видеть, что разработанный расчетный комплекс позволяет корректно вычислять гидродинамические характеристики крыльев даже предельно малой толщины.
1. Гурьев Ю. ВНикущенко Д. В. Опыт и перспективы использования компьютерных технологий для решения задач управляемости подводных объектов // Материалы науч.-техн. конференции «XII Макеевские чтения». СПб., 2005. С. 21—26.
2. Лифанов И. К. Метод сингулярных интегральных уравнений и численный эксперимент. М.: ТОО «Янус», 1995. 520 с.
3. Дворак А. В. Невырожденность матрицы метода дискретных вихрей в задачах пространственного обтекания // Труды ВВИА. 1986. № 1313. С. 441-453.
4. Лифанов И. К., Михайлов А. А. К расчету безотрывного и отрывного обтекания тел // Труды ВВИА. 1986. № 1313. С. 137-145.
5. Белоцерковский С. М., Ништ М. И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. М.: Наука, 1978. 352 с.
6. Гурьев Ю. В., Калинин О. С., Красиков В. И. и др. Нелинейная вихревая модель для определения гидродинамических характеристик кормового оперения ПЛ // Проблемы эксплуатации вооружения, военной техники и подготовки инженерных кадров ВМФ: сб. статей. СПб.: Изд-во ВМИИ, 2004. С. 54-56.
7. Никущенко Д. В., Калинин О. С., Красиков В. И., Рогожина Е. А. Определение гидродинамических характеристик изолированного оперения подводного аппарата на основе метода дискретных вихрей // 6-я Междунар. конференция по морским интеллектуальным технологиям «Моринтех-2005»: Труды конференции. СПб.: 2005. С. 75-76.
8. Никущенко Д. В., Рогожина Е. А. Влияние рулей на гидродинамические характеристики оперения подводных аппаратов // Науч.-техн. конференция «Крыловские чтения». 2006. С. 67.
9. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Лаборатория Базовых Знаний, 2001. 632 с.
10. Справочник по теории корабля: в 3 т. / под ред. Я. И. Войткунского. Л.: Судостроение, 1985. Т. 1: Гидромеханика. Сопротивление движению судов. Судовые движители. 768 с.
11. Белоцерковский С. М., Скрипач Б. К. Аэродинамические производные летательного аппарата и крыла при дозвуковых скоростях. М.: Наука, 1975. 424 с.
Статья рекомендована к печати проф. Ю. М. Далем.
Статья принята к печати 5 марта 2009 г.
