CC BY
466
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Моделирование обтекания тонкой пластинки с использованием модифицированной схемы метода вихревых элементов # 09, сентябрь 2013 DOI: 10.7463/0913.0602362
Макарова М. Е., Марчевский И. К., Морева В. С.
УДК 517.958:532.51
Россия, МГТУ им. Н.Э. Баумана masha-mak@mail.ru iliamarchevsky@mail.ru morevavs@rambler.ru
Введение
Моделирование обтекания тонкой пластинки потоком вязкой несжимаемой среды — классическая задача Блазиуса [1] — является одной из тестовый задач, используемык для верификации численных методов и алгоритмов вычислительной аэрогидродинамики.
В настоящей работе приведено решение этой задачи методом вихревыгх элементов. Отличие от метода, описанного в [2], состоит в использовании модифицированной расчетной схемы, в основу которой положено решение интегрального уравнения Фредгольма 2-го рода с ограниченным ядром. Такой подход позволил получить более точное решение поставленной задачи, правильно смоделировав профиль не только продольной (как представлено в [2]), но и поперечной скорости.
1. Постановка задачи
Рассмотрим модельную задачу об обтекании вязкой несжимаемой средой неподвижной тонкой пластинки (рис. 1), установленной вдоль набегающего потока. Будем считать, что скорость течения такова, что пограничный слой является ламинарным.
Течение среды описывается уравнениями Навье — Стокса [3]
dux dux dux 1 dp (d2ux д 2u + u--+ = —— + v -—r- +
dt x dx y dy pdx \ dx2 dy2
duy + u duy + u duy _ _ldp + /d2Uy + д2Uy\ ^
dt x dx y dy pdy \ дх2 dy° '
dux duy
—- +--- = 0,
dx dy
Рис. 1. Пограничный слой на обтекаемой пластинке
где V — коэффициент кинематической вязкости; р — плотность среды, которая полагается постоянной; и(г,Ь) = (их,иу)т — скорость; р = р(т,1) — давление в точке г = (х; у)т.
Уравнения (1) дополняются граничными условиями прилипания на поверхности обтекаемой пластинки и условиями затухания возмущений при удалении от нее.
Получение решения этих уравнений в пограничном слое вблизи пластинки основано на оценках-гипотезах о порядках различных членов в уравнениях (1) и пренебрежении малыми членами [3]. Уравнения пограничного слоя имеют вид:
dux
dt
+ Ux
dux dx
+ Uy
dux
d2ux 1 dp
dy dy2 p dux
dy
dux duy
(2)
0.
дх ду
Эти же уравнения можно также получить, исключая все малые слагаемые при предельном переходе в уравнениях (1) при Яе = и0¡/V ^ то, где Яе — число Рейнольдса [3], I — длина пластинки, и0 — скорость набегающего потока (рис. 1).
2. Решение Блазиуса
Будем считать пластинку бесконечно тонкой, а поток среды стационарным, u0 = const. Систему координат введем как указано на рис. 1: ось Ox направим вдоль пластинки, ось Oy — перпендикулярно к ней, начало координат поместим в переднюю кромку пластинки, тогда скорость установившего течения u = (ux; uy)T будет удовлетворять уравнениям [1, 3]:
du
ux
dx
x du + u
д2 ux
y
v-
dux +du
dy dy 0.
(3)
дх ду
Точное решение задачи (3) можно найти лишь для обтекания бесконечной пластинки (у = 0, х ^ 0), тогда граничные условия в приближенной постановке, предложенной Блазиусом, имеют следующий вид [1]:
ux = uy = 0 при y = 0, x ^ 0; ux ^ u0 при y ^ <х.
2
Решение задачи (3) с данными граничными условиями будет непригодно вблизи передней и задней кромок пластинки малой, но конечной толщины [1]. Если ввести функцию тока
дф дф ф(х,у) и* = ду, и = - дХ'
то уравнения (3) и указанные выше граничные условия можно привести к виду [1]:
дф д2ф дф д2ф д3ф
ду дхду дх ду2 ду3
дф дф
— = — = 0, ф = 0 при у = 0, х ^ 0, дх ду
(4)
дф ду
и° при у ^ то.
Уравнение (4) допускает автомодельное решение. В результате замены
п = у\—, / (п) = .-
у ух ^ихи°
оно приводится к обыкновенному дифференциальному уравнению
¿3/ , / 622 /
¿п3 + 2 ¿п
0
(5)
с краевыми условиями
/I = /
/ !п=о ¿п
п=о
0 ¿/ ' ¿п
(6)
п=<х
Краевую задачу (5)-(6) можно свести к задаче Коши относительно функции Г(С), где /(п)
а1/3р (а1/3п):
63Г Г ¿2Г + - — = 0,
¿С3 2 ¿С2
(7)
Г| = ¿Г ?=° ¿С
5=о
¿2С
0.
5=о
Эта задача с высокой точностью решается любым численным методом, при этом можно показать [3], что а & 0,332057.
Сделав обратные преобразования, возможно определить искомое стационарное поле скоростей в пограничном слое при обтекании тонкой пластинки.
3. Метод вихревых элементов
Численное решение задачи осуществляется с помощью метода вихревых элементов [2,4]. Достаточно тонкая плоская пластинка помещается в набегающий поток, который имеет
—>
2
1
скорость uo. По истечении некоторого времени течение устанавливается и определяется распределение скоростей в сечениях пограничного слоя на пластинке.
Скорость среды находится как суперпозиция полей скоростей отдельных вихревых элементов (ВЭ) и набегающего потока. ВЭ генерируются на поверхности обтекаемой пластины на каждом шаге расчета; их интенсивности таковы, что обеспечивается выполнение граничного условия на профиле [2, 4]. Таким образом, от правильности определения интенсивно-стей ВЭ зависит точность решения задачи в целом.
Для простоты выкладок течение считаем трехмерным с одинаковыми характеристиками во всех плоскостях z = const. Тогда u = (ux,uy, 0)T, интенсивность вихревого слоя на профиле в точке границы профиля дС с радиус-вектором ro = (xo,yo, 0)T также будем считать векторной величиной y(r0) = (0, 0,Y(ro))T. При этом скорость среды в точке с радиус-вектором r = (x,y, 0)T вычисляется по формуле
u(r) = Uo + Uwake(r) + / Y^) X (Г dlro , Г E R2 \ дС.
J 2n|r - ro|2
dC
Здесь uwake(r) — скорость среды, индуцируемая вихревым следом — ВЭ с известными положениями и интенсивностями, сгенерированными на предыдущих шагах расчета.
Предельное значение скорости потока u- (r0) на границе дС со стороны профиля равно [4]
u-(ro) = u(ro) - (X n(ro^ , ro E дС,
где n(ro) — единичная внешняя нормаль к профилю.
Для нахождения неизвестной плотности интенсивности вихревого слоя y(ro) в МВЭ используется условие равенства нулю вектора скорости u-(ro) на границе профиля, которое обеспечивается равенством нулю либо нормальной («классический» метод) [2, 4], либо касательной компоненты вектора u-(ro) («модифицированный» метод) [5, 6, 7]; с математической точки зрения эти подходы эквивалентны.
В «модифицированном» МВЭ решение задачи сводится к решению интегрального уравнения для y(ro):
Г [к X (r - ro)] -т(r)y(ro) diro - iM = -T(r) • (u^ + uwake), r E дС. (8) J 2n r—rJ 2
dC
Орт касательной т(г) к профилю выбирают таким образом, что п(г) х т(г) = к, где к — орт оси Ог.
Уравнение (8) — интегральное уравнение Фредгольма 2-го рода [4, 5], его ядро ограничено величиной к/4п, где к — кривизна границы профиля в соответствующей точке. Отметим, что «классический» МВЭ приводит к необходимости решения сингулярного интегрального уравнения с ядром Гильберта. В ограниченности ядра интегрального уравнения (8) и состоит основное преимущество «модифицированного» подхода над «классическим». В
связи с этим при его численном решении уравнения не накладывается жестких ограничений на вид расчетной схемы, она может быть достаточно произвольной [4, 5, 6, 7].
В расчете границу профиля разбивают на панели — отрезки ломаной, аппроксимирующей профиль; интенсивность вихревого слоя полагается постоянной на каждой панели и равной ее среднему значению 7г на этой панели. Тогда интегральное уравнение (8) может быть приближенно заменено системой из N линейных алгебраических уравнений, переход к которым подробно рассматривается в работе [7]
К
где Шгг0, М3Г0 — дифференциалы длин панелей с номерами г и ] соответственно, по которым производится интегрирование.
При численном решении системы (9) величины
и интегралы в правой части вычисляются аналитически по точным формулам [7]. Отметим, что в отличие от «классического» МВЭ выполнение граничного условия обеспечивается не в отдельных контрольных точках, а в среднем на каждой панели [5, 6, 7]. В связи с указанными особенностями использование модифицированного метода позволяет существенно повысить точность определения интенсивности вихревого слоя.
Отметим, что решение интегрального уравнения (9) не является единственным, поэтому необходимо ввести дополнительное условие
j) y(v)dlr = const,
dC
задающее циркуляцию скорости вдоль профиля, а получаемая система регуляризируется введением дополнительной переменной R [4].
Численное моделирование обтекания пластинки длиной I = 1,0, имеющей относительную толщину 2 %, передняя и задняя кромки которой имеют вид полуокружностей [2] осуществлялось средствами программного комплекса РОЬЛРЛ [8]. Профиль аппроксимировался ломаной, состоящей из 1276 панелей, безразмерная величина шага расчета по времени составляла Д£ = 1,5 • 10-4 при безразмерной скорости набегающего потока и0 = 1,0. Количество
Ki Kj
4. Результаты вычислительного эксперимента
ВЭ в расчетах превышало 100 000, поэтому активно использовались известные подходы к ускорению вычислений: параллельные алгоритмы и «быстрый метод» [8, 9, 10]. Задача решалась в нестационарной постановке, вязкость среды соответствовала значению Яе = 103. Скорость набегающего потока за первые 15 шагов расчета увеличивалась постепенно от нуля до значения и0, затем оставалась постоянной. Вид близкого к установившемуся вихревого следа вблизи пластинки представлен на рис. 2.
Рис. 2. Вид вихревого следа вблизи пластинки
В сечениях пластинки х = 0,35, х = 0,50, х = 0,65 измерялись величины скоростей. Для получения значения безразмерной скорости в расчете вместо значения и0 необходимо брать величину
тах и = щ (х) > и0
в соответствующем сечении, что связано с вытесняющим действием пластинки [2].
На рис. 3 в безразмерных величинах представлено распределение продольной компоненты скорости в сечении х = 0,35.
Рис. 3. График распределения приведенной продольной компоненты
скорости в сечении х = 0,35 пограничного слоя на пластинке (точки — результаты расчета, сплошная линия — решение Блазиуса)
Хорошее совпадение результатов расчета со значениями, полученными аналитически, наблюдалось во всех рассматриваемых сечениях пластинки (х = 0,35, х = 0,50, х = 0,65) во все моменты времени после того, как течение вышло на стационарный режим.
Поперечная компонента скорости иу по величине существенно меньше и не может быть выражена в безразмерных переменных. Однако известно, что [1]
иу
ио
дф дх
2'
^ (/' - /).
Поэтому можно провести анализ величины
и
x
и0 у vu0
Uy
Uo
x Re
l
и построить для нее график в безразмерных переменных Блазиуса (рис. 4). Хорошее согласие результатов наблюдается также и в других сечениях x = const.
Рис. 4. График распределения поперечной компоненты скорости в сечении х = 0,35 пограничного слоя на пластинке (точки — результаты расчета, сплошная линия — решение Блазиуса)
Отметим, что профиль продольной скорости, показанный на рис. 3, соответствует мгновенному распределению, тогда как профиль поперечной скорости существенно меняется от шага к шагу; представленный на рис. 4 график получен путем осреднения по большому числу шагов (от 100 шагов в сечении х = 0,35 до 300 шагов в сечении х = 0,65).
Расчеты показывают, что при использовании «классического» МВЭ, приводящему к необходимости решения сингулярного интегрального уравнения на профиле, продольная компонента скорости в установившемся течении также хорошо согласуется с аналитическим решением, а для получения правильного распределения поперечной скорости требуется значительно увеличивать число ВЭ в расчетной схеме и/или проводить осреднение результатов по значительно большему временному промежутку.
Заключение
В работе рассмотрена модельная задача Блазиуса о моделировании методом вихревых элементов обтекания тонкой пластинки конечной длины вязкой несжимаемой средой. Интенсивности вихревых элементов, моделирующих пограничный слой, определялись «модифицированным» методом МВЭ с касательными компонентами скорости, реализованным в программном комплексе РОЬЛЯЛ. Проведенное сравнение полученных результатов с известным аналитическим решением Блазиуса показало высокую точность предложенной «модифицированной» расчетной схемы МВЭ.
Работа выполнена при финансовой поддержке гранта Президента РФ для государственной поддержки молодых российских ученых — кандидатов наук (проект МК-6482.2012.8) и гранта Президента РФ для государственной поддержки ведущих научных школ (проект НШ-255.2012.8).
Список литературы
1. Лойцянский Л.Г. Ламинарный пограничный слой. М.: ГИФМЛ, 1962. 480 с.
2. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. М.: Изд-во МГУ, 2006. 184 с.
3. Седов Л.И. Механика сплошной среды. СПб.: Лань, 2004. 560 с.
4. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М.: Янус, 1995. 520 с.
5. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H. Accuracy Considerations for Implementing Velocity Boundary Conditions in Vorticity Formulations // SANDIA Report SAND96-0583 UC-700. March, 1996. 50 p.
6. Макарова М.Е. Расчет стационарного безотрывного обтекания профиля потоком идеальной несжимаемой среды // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2011. Спец. выпуск «Прикладная математика». C. 124-133.
7. Морева В.С. Вычисление вихревого влияния в модифицированной схеме метода вихревых элементов // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2012. Спец. выпуск №2 «Математическое моделирование в технике». C. 137-144.
8. Марчевский И.К., Морева В.С. Численное моделирование в задачах аэрогидродинамики с использованием метода вихревых элементов // CAD/CAM/CAE Observer. 2012. №2. С. 84-91.
9. Морева В.С. Способы ускорения вычислений при решении плоских задач аэродинамики методом вихревых элементов // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2011. Спец. выпуск «Прикладная математика». C. 83-95.
10. Лукин В.В., Марчевский И.К., Морева В.С., Попов А.Ю., Шаповалов К.Л., Щеглов Г.А. Учебно-экспериментальный вычислительный кластер. Ч. 2. Примеры решения задач // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. 2012. №4. С. 82-102.
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Flow simulation around a thin plate using a modified numerical scheme of the vortex element method
# 09, September 2013
DOI: 10.7463/0913.0602362
Makarova M. Ye., Marchevsky I. K., Moreva V. S.
Bauman Moscow State Technical University 105005, Moscow, Russian Federation masha-mak@mail.ru iliamarchevsky@mail.ru morevavs@rambler.ru
The Blasius benchmark problem dedicated to numerical simulation of a flow around a thin plate is considered in this paper. The laminar boundary layer was simulated using a modified numerical scheme of the vortex element method with tangent velocity components. This numerical scheme was implemented in POLARA software. Distributions of longitudinal and transversal velocity components over the boundary layers cross-section were obtained. Use of the modified scheme allowed to improve accuracy of determining the intensity of a vortex layer on the plate significantly at every step of the calculation. The obtained results appeared to be in good agreement with the analytic solution given by Blasius.
References
1. Loytsyanskiy L.G. Laminarnyypogranichnyy sloy [Laminar boundary layer]. Moscow, GIFML, 1962. 480 p.
2. Andronov P.R., Guverniuk S.V., Dynnikova G.Ia. Vikhrevye metody rascheta nestatsionarnykh gidrodinamicheskikh nagruzok [Vortex methods of calculation of unsteady hydrodynamic loads]. Moscow, MSU Publ., 2006. 184 p.
3. Sedov L.I. Mekhanika sploshnoy sredy [Mechanics of continuum]. St. Petersburg, Lan', 2004. 560 p.
4. Lifanov I.K. Metod singulyarnykh integral'nykh uravneniy i chislennyy eksperiment [Method of singular integral equations and numerical experiment]. Moscow, Yanus, 1995. 520 p.
5. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H. Accuracy Considerations for Implementing Velocity Boundary Conditions in Vorticity Formulations. SANDIA Report SAND96-0583 UC-700. Sandia National Laboratories, Albuquerque, N.M., March 1996. 50 p.
6. Makarova M.E. Raschet statsionarnogo bezotryvnogo obtekaniya profilya potokom ideal'noy neszhimaemoy sredy [Calculation of stationary attached flow past airfoil by stream of ideal incompressible medium]. VestnikMGTUim. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2011, spets. vyp. "Prikladnaya matematika" [spec. iss. "Applied mathematics"], pp. 63-74.
7. Moreva V.S. Vychislenie vikhrevogo vliyaniya v modifitsirovannoy skheme metoda vikhrevykh elementov [Calculation of vortex effect in the modified numerical scheme of vortex-element methods]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2011, spets. vyp. "Matematicheskoe modelirovanie v tekhnike" [spec. iss. "Mathematical modeling in engineering"], pp. 137-144.
8. Marchevskiy I.K., Moreva V.S. Chislennoe modelirovanie v zadachakh aerogidrodinamiki s ispol'zovaniem metoda vikhrevykh elementov. CAD/CAM/CAE Observer, 2012, no. 2, pp. 84-91.
9. Moreva V.S. Sposoby uskoreniya vychisleniy pri reshenii ploskikh zadach aerodinamiki metodom vikhrevykh elementov [Ways for calculation speed-up in solving 2D aerodynamics problems by vortex element method]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2011, spets. vyp. "Prikladnaya matematika" [spec. iss. "Applied mathematics"], pp. 83-95.
10. Lukin V.V., Marchevskiy I.K., Moreva V.S., Popov A.Yu., Shapovalov K.L., Shcheglov G.A. Uchebno-eksperimental'nyy vychislitel'nyy klaster. Ch. 2. Primery resheniya zadach [Computing cluster for training and experiments. Part 2. Examples of solving problems]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2012, no. 4, pp. 82-102.
