CC BY
6
УДК 531.6.011.32: 532.582.4 БОТ: 10.18698/2309-3684-2021-4-80102
Методика организации итерационного процесса при моделировании дозвукового отрывного обтекания удлиненных тел с применением эквивалентной поверхности и кубических сплайнов
© В.Н. Тимофеев МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Рассмотрены вопросы организации итерационного процесса построения поверхности эквивалентного тела в задаче математического моделирования дозвукового отрывного обтекания удлиненных тел с частичной реализацией концепции вязко-невязкого взаимодействия. Использовалась схема течения с полубесконечной эквивалентной поверхностью. Численное моделирование проведено по алгоритмам методики с применением метода дискретных вихрей. Для большей универсальности при построении расчетной сетки на поверхности эквивалентного тела использовалась аппроксимация сглаживающими кубическими сплайнами. Представлены данные о влиянии формы хвостового участка эквивалентной поверхности на распределение скорости и давления при осесимметричном обтекании тел с донным срезом.
Ключевые слова: математическое моделирование, дозвуковое отрывное обтекание, концепция вязко-невязкого взаимодействия, эквивалентная поверхность, сглаживающий кубический сплайн
Введение. Для обширной сферы практических приложений важность, актуальность и сложность проведения математического моделирования дозвукового отрывного обтекания пространственных тел породили множество подходов, методов и методик [1-11] к решению данной задачи. При исследовании обтекания с отрывом потока на удлиненных телах с донным срезом в работе [12] была представлена методика, построенная на основе частичной реализации концепции вязко-невязкого взаимодействия. Основными целями, при создании данной методики являлись получение распределения скоростей и давлений по поверхности обтекаемого тела и снижение трудоемкости вычислений для нахождения таких распределений. Рассматривались варианты реализации методики для двух схем течения: с эквивалентным телом конечной длины [12] и с полубесконечным эквивалентным телом [13]. В обеих схемах течения поверхность эквивалентного тела называется эквивалентной и состоит из носовой ^ и хвостовой Е2 частей. Геометрия носовой части эквивалентной поверхности определяется формой обтекаемого тела, а хвостовая часть Е2 моделирует влияние следа, который развивается за обтекаемым телом в реальном вязком течении. В случае рассмотрения
режимов течения, для которых область отрыва потока расположена вблизи контура донного среза, предполагается, что линия отрыва совпадает с этим контуром. При исследовании удлиненных тел считается, что толщина вытеснения пограничного слоя пренебрежимо мала и носовая часть ^ эквивалентной поверхности совпадает с
поверхностью обтекаемого тела. Формирование поверхности Е2 хвостовой части эквивалентного тела осуществляется в ходе проведения итерационного процесса. В работе [14] излагаются общие сведения об организации данного процесса. В настоящей статье порядок проведения итераций излагается в более развернутом форме. Частично используются терминология и обозначения работы [14].
Начальная стадия итерационного процесса. Рассматривается дозвуковое обтекание осесимметричных тел. В качестве начальных данных задаются конфигурация поверхности обтекаемого тела и параметры набегающего потока, в том числе, его вектор скорости ¿>ш . Используются безразмерные длины Ь и координаты х, у, г прямоугольной декартовой системы координат Оху (рис. 1), связанной с обтекаемым телом. Безразмерные параметры получаются отнесением значений соответствующих размерных величин к диаметру миделевого сечения Б* . Поэтому безразмерный диаметр миделевого сечения оказывается равным единице: Бт = 1.
Рис. 1. Поверхность эквивалентного тела: носовая ^ и хвостовая Е2 части эквивалентной поверхности; Ь. — линия их стыковки; Еп и Е12 — головной и цилиндрический участки
носовой части с длинами Ь и Ь ; Е21 — участок хвостовой части длиной Ь5 и переменным диаметром Б (х); ^22 — полубесконечный участок с постоянным диаметром Б ; Ь — линия соединения участков Е21 и Ъгг
Исследуемые тела состоят из оживальной головной части с длиной Д , присоединенной к цилиндрической части, имеющей длину
Д и плоский донный срез. В соответствии с такой геометрией носовая часть ^ эквивалентной поверхности представляет собой объединение двух участков оживальной £п и цилиндрической £12 формы. Конфигурация поверхности обтекаемого тела и поверхностей обоих участков £п и £12 в полной мере определяется безразмерными
длинами Д и £ [15]. Сумма Ь + Д обозначается через Д и
равняется полной длине обтекаемого тела и носовой части эквивалентной поверхности. Исследуется обтекание удлиненных тел, у которых длина цилиндрической части более чем в два раза превышает длину
головной части (Д > 2Д ) .
Применяется схема течения с полубесконечным эквивалентным телом. Носовая часть ^ эквивалентной поверхности соединяется с ее хвостовой частью Е2 по линии стыковки Д . Согласно данным работы
[3] поперечный размер полубесконечного тела вытеснения, моделирующего влияние следа, является переменным и принимает на бесконечности некоторое фиксированное значение. В используемой методике это значение обозначается через Д и считается, что оно
достигается при конечном удалении Д от донного среза. Таким образом, на удалениях от донного среза, не превышающих Д, хвостовая часть эквивалентного тела имеет переменный диаметр Д (х), а при больших удалениях — постоянный диаметр Д . Поэтому хвостовая часть Е2 эквивалентной поверхности образуется соединением по линии Д участка Е21 конечной длины Д с переменным диаметром Д (х) и полубесконечного участка Е22, который имеет постоянный диаметр Д. Параметры Д и Д являются важными геометрическими характеристиками участка Е21. Для краткости будем называть их длиной и минимальным диаметром следа. При исследовании осесимметричного обтекания с нулевым углом атаки для полного определения конфигурации поверхности Е2
хвостовой части эквивалентного тела на участке Е21 достаточно задать зависимость Д (х) диаметра Д поперечного сечения от продольной координаты х . Проведение математического моделирования осуществляется для трех вариантов зависимости Д (х).
Задание длин Ь и Ьс участков носовой части эквивалентной по-
Т(0)
верхности , начальных значений длины ц' и минимального диаметра следа, а также выбор одного из вариантов Б(0) (х) зависимости Б (х) позволяет полностью определить геометрическую форму эквивалентной поверхности в начальном приближении при осесимметрич-ном обтекании исследуемого тела. Затем с помощью аппроксимации сглаживающими кубическими сплайнами [16] производится построение расчетной сетки Е(0) на обоих участках Еп и Е12 носовой части эквивалентной поверхности. Эта сетка, как и значения длин Ь и Ц
указанных участков, не изменяются при дальнейших вычислениях. Аналогичным образом с применением сплайн аппроксимации, базиру-
г(0)
ясь на заданных начальных значениях длины Ь и минимального диаметра Б^ следа и на выбранном варианте зависимости Б20 (х), строится начальное (нулевое) приближение расчетной сетки Е(20) на участках Е21 и Е22 хвостовой части эквивалентной поверхности.
Объединением двух построенных сеток получается начальное (нулевое) приближение расчетной сетки для всей эквивалентной поверхности Е(0) = Е(0) и Е20).
Первый этап итерационного процесса. В соответствии с основным положением концепции вязко-невязкого взаимодействия скорости и давления на поверхности обтекаемого тела определяются по результатам расчета невязкого обтекания эквивалентного тела. Так как форма эквивалентной поверхности первоначально задается только начальным приближением Е(0), в применяемой методике для окончательного формирования поверхности эквивалентного тела организуется итерационный процесс. На каждой итерации будем рассматривать несколько этапов. В свою очередь могут выделяться составные части этапа, которые будем называть стадиями.
На первом этапе у - ой итерации производится формирование
поверхности Е^' такого эквивалентного тела, для которого на втором этапе проводится расчет невязкого обтекания. Будем называть Е^ эквивалентной поверхностью начального приближения для у - ой итерации. На первой итерации (у = 1) в качестве поверхности Е(1)
берется поверхность Е(0) = Е(0) и Е(0), описание построения которой было проведено выше. Естественно, что при этом длина и минимальный диаметр следа для эквивалентной поверхности Е начального
приближения в первой итерации задаются следующим образом: Д!)= ДО) и Д) = Д(°) .
В последующих итерациях (] > 2) поверхность ) формируется по результатам расчетов на предыдущей итерации. А именно, по значениям длины Д -1) и минимального диаметра Д7-1) следа, по виду
зависимости 1 (х), используемой для построения эквивалентной
поверхности 1 начального приближения в (7 -1) - ой итерации, и по результатам расчета невязкого обтекания эквивалентного тела с поверхностью 1 производится вычисление величин Д), Д7 и определение вида зависимости Д7 (х) для 7 - ой итерации. Затем с применением сглаживающих кубических сплайнов производится
^ 7)
построение расчетной сетки для эквивалентной поверхности 10 начального приближения в 7 - ой итерации. Более подробное изложение этих действий, определяющих геометрию эквивалентной поверхности начального приближения для 7 - ой итерации, дается ниже при описании второго и третьего этапов данной итерации.
Здесь и ниже, во избежание громоздкости верхними индексами (7), указывающими на принадлежность к действиям, производимым в ходе 7 - ой итерации, снабжены лишь некоторые из рассматриваемых объектов и величин, например, эквивалентная поверхность
1(7). В случаях, когда ясно, что действия осуществляются на 7 - ой
итерации, остальные объекты и величины индексами (7)
не снабжаются.
Второй этап итерационного процесса. На втором этапе каждой
итерации производится численное моделирование невязкого
обтекания тела, поверхность которого представляет собой эквивалент-
ную поверхность начального приближения для 7 - ой итерации. При проведении моделирования предполагается, что газ, обтекающий исследуемое тело, является идеальным и невесомым. Для описания динамики такой газовой среды используются уравнения неразрывности и движения в форме Эйлера. Требуется, чтобы решение системы, составленной из этих уравнений, удовлетворяло граничным условиям непротекания поверхности обтекаемого тела и затухания возмущений на бесконечности. Рассматривается установившееся движение газовой среды с умеренными дозвуковыми скоростями набегающего потока. Это позволяет считать среду несжимаемой и однородной.
Постулируется, что в области пространства, расположенной во внешности поверхности эквивалентного тела, течение потенциально и существует потенциал р-р(х,у,2) вектора скорости возмущений.
Тогда скорость газового потока в произвольной точке М0(х0, у0, ^)
указанной области находится через градиент потенциала вектора скорости возмущений:
и(Мй) = и„+Ч<р(Мй), (1)
— <9 — д — д _ _
где У = е,--\-е2--\-еъ--оператор Гамильтона [17]; ех,е2,еъ —
Эх0 Эх0 Эх0
базисные орты выбранной системы координат Охуг .
При сформулированных выше условиях потенциал р вектора скорости возмущений является решением внешней задачи Неймана для уравнения Лапласа. Граничные условия в задаче Неймана получаются из граничных условий непротекания поверхности обтекаемого тела и затухания возмущений на бесконечности, соблюдение которых требуется для системы из уравнений неразрывности и движения газовой среды. Решение внешней задачи Неймана представляется в виде потенциала двойного слоя, что обеспечивает потенциалу вектора скорости возмущений удовлетворение уравнения Лапласа и условий затухания возмущений на бесконечности.
Удовлетворение оставшегося граничного условия непротекания производится следующим образом. Записывается соотношение для нахождения градиента потенциала вектора скорости возмущений:
где Е — эквивалентная поверхность Е^ начального приближения для у — ой итерации; йо — площадь текущего элемента й Е эквивалентной поверхности Е; 7 — вектор, соединяющий точку М0 вычисления параметров с точкой М , которая находится на элементе поверхности ¿/Е; г — модуль вектора 7; я(М) — орт внешней
нормали к поверхности Е, построенный в точке М; g(M) —
поверхностная плотность потенциала двойного слоя.
Знание градиента потенциала вектора скорости возмущений в случае, когда точка М находится во внешности поверхности эквивалентного тела, позволяет вычислить вектор скорости газового потока по соотношению (1). Если же точка М0 принадлежит
поверхности эквивалентного тела , на которой задается поверхностная плотность g(М), то градиент М0) потенциала двойного слоя испытывает в этой точке разрыв [18]
8 О7Ф) (мо) = (м0) т, (м0)+§'Т2 (м0) (м0),
где (М) и (М) — частные производные поверхностной плотности, вычисляемые по направлениям взаимно ортогональных ортов тх (М0) и т2 (М{)) в плоскости, проходящей через точку М0 и касательной к поверхности 2.
Наличие разрыва с>(У$>)(М0) приводит к тому, что предельное значение градиента потенциала двойного слоя, вычисляемое в точке М0 на поверхности эквивалентного тела, оказывается равным сумме
Последнее выражение используется при нахождении вектора скорости потока газа в точках поверхности эквивалентного тела:
й (М „) = В.,V|| Г/Н(М) х)с1о+ Х-сЪ(М,,), (3)
где Зи(М0) — вектор разрыва скорости потока, равный вектору
с>(У<^)(М0) разрыва градиента потенциала вектора скорости возмущений.
Если форма и геометрические параметры поверхности 1 известны, то при вычислении интегралов в соотношениях (2) и (3) могут быть определены все величины за исключением поверхностной плотности потенциала двойного слоя. В рассматриваемой методике нахождение поверхностной плотности потенциала двойного слоя, а следовательно, и градиента потенциала вектора скорости возмущений и скорости газового потока, выполняется с помощью численного моделирования [19]. Для этого на одной из стадий первого этапа 7 - ой итерации осуществляется построение расчетной сетки. Участки
, Е22 и Е21 эквивалентной поверхности ) начального приближения для у - ой итерации с применением сглаживающих кубических сплайнов аппроксимируется набором панелей 1, к=1,...,N .
Каждая из этих панелей имеет форму либо четырехугольника, либо треугольника. Как было указано выше, геометрия участков 1И и
Е22 носовой части эквивалентной поверхности не изменяется в ходе
проведении итерационного процесса. Поэтому геометрия панелей, которые располагаются на этих участках также остается неизменной:
Е ^ ^ =Е. Поверхность полубесконечного участка Е22 эквивалентного тела заменяется совокупностью полубесконечных панелей Е^', к=Ж+1,..., N+N. Расчетная сетка получается путем объединения всех аппроксимирующих панелей Е[3^, к=1,..., N+N. Эта сетка образует
Е ^ эквивалентную поверхность начального приближения для 3 — ой итерации.
При численном моделировании плотность потенциала двойного слоя на любой панели из сетки считается постоянной величиной: g(M) = gi:, к = 1,...,N+N. Использование расчетной сетки позволяет
преобразовать поверхностные интегралы в правых частях равенств (2) и (3) в суммы поверхностных интегралов, вычисляемых по всем панелям Е3, к=1,..., N+N.
Далее учитывается, что градиент потенциала двойного слоя,
„ -А(3)
размещенного на каждой из панелей Е^ ' с постоянной поверхностной плотностью ^, оказывается равным вектору скорости, индуцированной вихревой линией Ьк, расположенной на границе сЕ^ к — ой панели, при условии, что указанная вихревая линия имеет циркуляцию
Гк =—gk, k = 1,...,N+Nr. (4)
Поэтому на границах ^ всех панелей Е3, к=1,..., N+N размещаются вихревые многоугольники ^. На участках Еи, Е22 и Е поверхности эквивалентного тела размещаются вихревые четырехугольники или треугольники, а на полубесконечном участке Е22 — П—образные вихревые линии. Векторы скорости,
индуцируемые в данной точке потока вихревым линиями, находятся с использованием формулы Био-Савара. Для большей компактности записи расчетных соотношений при вычислениях по формуле
Био-Савара используются векторы называемые векторами
функции скорости, индуцированной в точке М0 к — ой вихревой линией, и определяемые по формуле:
—, Л 1 гЛхг
^ Мо — — (5)
4^ г
ьк
Вычисление векторов функций скорости м^(М0), к /V + ,
производится по алгоритмам [20], адаптированным для проведения численного моделирования. Предполагается, что вихревые многоугольники конечного размера и П -образные вихри, расположенные по разные стороны от линии Ьи соединения участков Е21 и Е22, имеют равные циркуляции. Это дает возможность сократить число вихревых образований путем объединения указанных пар вихревых линии с нахождением преобразованных векторов функций скорости
щ. (М0), к = 1,..., /V . В конечном итоге получается, что, если точка М0
располагается в поле течения, например, во внешности поверхности эквивалентного тела, то вектор скорости газового потока определяется по следующему соотношению:
N _.
й{М0) = йх+^Гкм?1{М0), (6)
к=1
Если же точка М0 находится на поверхности эквивалентного тела, то
формула для вектора скорости газового потока приобретает следующий вид:
N -> 1 _,
й(М„) = й.,+%Г1с м>1 (М0) + - би (М0). (7)
к=1 2
Представление вектора скорости равенствами (6) и (7) позволяет использовать алгоритмы метода дискретных вихрей. В соответствии с данными алгоритмами граничное условие непротекания должно выполняться на множестве контрольных точек С, ^ = !>•••> N
расположенных на панелях Е. Это позволяет найти неизвестные циркуляции Гк, k=1,..., N как решение следующей системы линейных алгебраических уравнений:
N
Еа*кГк = К, У = 1.....N, (8)
к=\
где аук=щ(Су)-п(Су)-, Ъу = -п(Су); у = 1,...,М; к = \,...,М; Я(СУ) — орт внешней нормали к поверхности панели построенный в контрольной точке С„.
Решение системы (8) производится с использованием регуляризи-рующей функции [21]. После этого векторы скорости потока в точках поля течения находятся по равенству (6). Равенства (4) позволяют найти на всех панелях величины поверхностной плотности ^, к=1,...,N+N . Поэтому распределение скорости потока газа на
поверхности эквивалентного тела получается с применением соотношения (7) для вычисления векторов скорости потока в контрольных точках:
N -. 1 _.
= + У = 1,..,М. (9)
к-1 2
Затем в контрольных точках определяются значения модуля вектора скорости газового потока и = |и(Су)| и безразмерные
скорости и=и/игде и., =\о.,\ — модуль вектора скорости
набегающего потока.
Для определения статического давления в точках поля течения и поверхности эквивалентного тела используется интеграл Бернулли
р - р» + 1 лД^-о2), (10)
где р и о — статическое давление и модуль вектора скорости газовой среды в данной точке; рш и л — статическое давление и плотность набегающего потока. Применение интеграла Бернулли (10) позволяет
2 (р — рш )
найти распределение коэффициента давления с - —-на поЛ Ооо
верхности эквивалентного тела по формуле
Ср - 1 — (о)2. (11)
Третий этап и порядок проведения итерационного процесса.
Расчет невязкого обтекания эквивалентного тела с поверхностью ) выполняется на втором этапе каждой ' - ой итерации с применением соотношений (1)—(11). На третьем этапе результаты этого расчета используются для коррекций значений длины и минимального диаметра следа, а также вида зависимости Д (х) диаметра поперечного сечения на участке поверхности Е21. Вычисление скорректированных значений указанных параметров Д;+1) и Д'+1), и получение уточненной зависимости Д('+1) (х) осуществляются следующим образом.
На завершающей стадии второго этапа текущей итерации значения коэффициента давления с находятся в точках пересечения нескольких продольных сечений эквивалентной поверхности и линии Lj стыковки носовой и хвостовой частей, совпадающей с контуром
донного среза обтекаемого тела. Эти значения получаются
интерполяций по величинам коэффициента давления, вычисленных в контрольных точках, лежащих на данном продольном сечении по разные стороны от указанных точек пересечения. Затем осреднением по контуру Lj величин коэффициента давления в рассмотренных
точках пересечения определяется расчетное значение коэффициента донного давления (с^ ) .
Далее расчетное значение коэффициента донного давления (с^ )(^) сравнивается с оценочным значением коэффициента донного
давления (ср!,) . Последнее получается по материалам работ [13,19] с
использованием формулы Хорнера, обобщающей экспериментальные данные об аэродинамическом коэффициенте донного сопротивления
удлиненных тел. Для исследуемых тел величина (с^ ) не изменяется
в ходе итерационного процесса, так как является функцией суммы безразмерных длин Д и Lс, а также, числа Рейнольдса Яе, которое находится по модулю скорости набегающего потока и полной длине обтекаемого тела £ь.
Разность расчетного и оценочного значений коэффициента донного давления определяет получаемое линейной интерполяцией промежуточное значение приращения для длины следа
= , У ° ' } А(с Т-(с ) . (12)
' (ЫГ-ыгт с е>
Здесь индексом ' обозначены соответствующие величины на 1 - ой итерации; индексом ( ' -1 ) при дополнительном условии о том, что 1 > 2 , снабжаются параметры на ( 1 -1) - ой итерации; для первой
итерации ' -1 - 0 и следует положить - 0 и (с^)( ' - 0 .
Окончательное значение приращения и уточненная величина длины следа определяются на первом этапе (1 +1 ) - ой итерации:
(13)
(14)
где Л( 1+1) — коэффициент релаксации.
Похожий алгоритм используется для определяется корректировки
значения Д1 минимального диаметра следа. Второй этап текущей
итерации завершается нахождением аэродинамического
коэффициента суммарного сопротивления обтекаемого тела С, равного сумме аэродинамических коэффициентов сопротивления давления С и трения и коэффициента донного сопротивления
{рх* ){'. Аэродинамические коэффициенты сопротивления давления и трения определяются интегрированием по носовой поверхности эквивалентного тела коэффициентов давления е и местного
коэффициента поверхностного трения е находится по формуле для турбулентного пограничного слоя на плоской пластине
е
= 0,455{^Яе) 258, где Яе — число Рейнольдса.
Для рассматриваемых в данной статье тел, у которых цилиндрический участок носовой части эквивалентной поверхности имеет плоский донный срез, значение коэффициента донного сопротивления
(е**)' равняется взятому с обратным знаком расчетному значению
коэффициента донного давления (е^. В соответствии с данными работы [12] величина минимального диаметра должна равняться
значению ^' = ^0,5е{.') . Поэтому предварительное значение приращения для минимального диаметра следа определяется как разность
(15)
Окончательное значение приращения и уточненная величина длины следа для первого этапа ( ' +1) - ой итерации определяются по следующим формулам:
(16)
= £>(;) +М)[1+1\ (17)
где — коэффициент релаксации . При проведении численного
моделирования часто используются следующие пары значений коэффициентов релаксации: = 0,6...1; = 0 и = 0; ^'+1)= 0,6...1.
Вычислительный эксперимент проводится для трех вариантов выбора функции Б 0' (х), каждый из которых определяет зависимость
диаметра поперечного сечения хвостовой части эквивалентного тела от продольной координаты х в начальном приближении. В первом варианте функция Б0) (х) получается путем ряда преобразований косинусоиды:
D = D
(0)
(x) = 0,5 (1 + Di°)) + 0,5 (1 - D<°)) cos *( x - L)
V
(18)
Здесь Д - Д + Д — полная длина обтекаемого тела; и Д0,1 —
VI1)
длина и минимальный диаметр следа эквивалентной поверхности Еу в начальном приближении. Для второго варианта выбирается функция
А Г , хЬ ЛЛ
)(0)
D = D
(0)
(x ) = 1 -(1 - Dl" >)
1 -
1-
f Т Л
x - L
V V
L(
0)
exp (-c(x - Lb))
(19)
JJ
при задании которой параметр b определяется по формуле b = Vк 2, где к — константа, значения которой выбираются в диапазоне (0,1...0,4); параметр c = -b/L1 . В соотношениях (18) и (19) значения
продольной координаты x изменяются от Lb до Lb + L(0). Третий вариант получается из второго отбрасыванием ряда точек, для которых продольная координата x принимает значения близкие к L + L(0). Естественно, что в этом случае длина L(0) определяется с учетом длины отрезка, содержащего отброшенные точки.
При проведении первой итерации зависимость D(1) (x) совпадает
с зависимостью D(0) (x). В каждой последующей итерации (при j > 2 ) на её первом этапе зависимость D(j) (x) диаметра поперечного сечения хвостовой части от продольной координаты x - ой итерации получается преобразованием зависимости D(j-1) (x), которая использовалась на (j -1) - ой итерации. Для этого с учетом величин
Lj) и Dä(j), , вычисляемых по соотношениям аналогичным (12)-(17) но с уменьшенным на единицу значением индекса, выполняются два последовательных преобразования:
D = D{]){x) = D(]-l)\ (:x-Lb)
L:
L\J
D = D[1) (x) =
(1 - D.(j))
x - L
(1 -Dij-1))V Lj J
D(]\x).
(20)
(21)
Итерационный процесс завершается после достижения заданной точности определения значений длины Д ^ и минимального диаметра
следа Д'. Конечные величины длины и минимального диаметра
следа обозначаются через Де) и Д(е).
Результаты вычислительного эксперимента. В качестве иллюстрации некоторых важных моментов, связанных с проведением итерационного процесса при математическом моделировании дозвукового отрывного обтекания пространственных тел по методике с частичной реализацией концепции вязко-невязкого взаимодействия на базе применения схемы течения с полубесконечным эквивалентным телом и расчетных соотношений (1)-(21), представляются результаты вычислительного эксперимента. Вычислительный эксперимент проводился для случая обтекания осесимметричного оживально-цилиндрического тела с нулевым углом атаки при числе Рейнольдса Яе = 5 -106. У исследуемого тела длины головной оживальной и цилиндрической частей равняются Ьг = 1 и = 2 .
Рис. 2. Окончательная форма эквивалентной поверхности при обтекании тела Ь = 1 и ¿г = 2 при а = О и Ке = 5 • 10'1 для первого варианта функции
Д0) (х): Де) = 2,53 и Д{е) = 0,291
Такая форма поверхности эквивалентного тела получается при длине следа, равной Де) = 2,53, и при минимальном диаметре следа, равном = 0,291. На рис. 2 и на последующих рисунках не показывается полубесконечный участок Е22 поверхности эквивалентного тела с постоянным значением диаметра Д и данные о расчетных параметрах на этом участке. Начальная и окончательная конфигурации поверхности эквивалентного тела для второго варианта задания
функции Д(х), который соответствует формуле (19), изображаются на рис. 3.
Начальная конфигурация получается при подстановке в зависимость (19) начальных значений длины следа = 5,600 и
минимального диаметра следа Д0) = 0,280 и последующего
построения расчетной сетки Е(0). У полученной после завершения итерационного процесса окончательной конфигурации эквивалентной поверхности конечное значение длины следа оказывается равным
, а минимального диаметра следа — = 0,291. Таким образом, в ходе проведения итераций расчетная длина следа Ь5 изменяется существенным образом. Расчетное значение минимального диаметра следа Б,, также изменяется, но не столь существенно как . На эквивалентных поверхностях, изображенных на рис. 3, дополнительно с помощью цветовых диаграмм представлено распределение модуля безразмерной скорости и газового потока, которое находится в контрольных точках вычислениями по формуле (9). Для третьего
варианта задания зависимости 0) (х) начальная и окончательная
конфигурации поверхности эквивалентного тела изображаются на рис. 4.
б
Рис. 3. Формы эквивалентной поверхности и распределение безразмерной скорости П при обтекании тела с!г=1 и!с=2 при а = 0 и Яе = 5• 106: а — начальная конфигурация для второго варианта функции Б(0)(х) с Д0) = 5,600 и Д(0)= 0,280 ;
б — окончательная конфигурация с ¿е) = 4,617 и Б(0) = 0,291
а
Рис. 4. Формы эквивалентной поверхности и распределение коэффициента давления с при обтекании тела с ¿„ = I и £ = 2 при а = 0 и Ке = 5 • 106: а — начальная конфигурация для третего варианта функции Д0) (х) с Д0)- 3,400 и Д0)- 0,277 ;
б — окончательная конфигурация с Де) - 3,197 и Д0) - 0,291
При начальных значениях длины следа - 3,400 и минимального диаметра следа Д0) - 0,277 для третьего варианта задания функции Д0' (х) конечные значения этих параметров оказываются
равными следующим величинам: -3,197 и Д0) -0,291. Это
свидетельствует о том, что, как и для второго варианта в ходе проведения итераций расчетная значения длины и минимального диаметра следа претерпевают заметные изменения. На рис. 4 также представлены цветовые диаграммы, показывающие распределение по поверхности эквивалентного тела коэффициента давления с , найденное в контрольных точках из интеграла Бернулли по формуле (11).
Так как при использовании данной методики одной из основных задач является получение распределений скорости и давления по поверхности обтекаемого тела, то в продольном сечении поверхности эквивалентного тела находятся зависимости безразмерной скорости й = о(х) и коэффициента давления с = с (х) от координаты х.
Графики указанных зависимостей представлены на рис. 5а и 5 б, соответственно, для первого (кривые 1) и второго (кривые 2) вариантов выбора функции Б(0) (х).
0,8
0,4
0,0
Г 2
1
0,8
0,4
0,0
-0,4
б
Рис. 5. Зависимости параметров потока газа в продольном сечении поверхности эквивалентного тела от координаты х при обтекании тела с Ь =1 и Ьс =2 при а = 0 и Ке=5-106:
£> = 2,53 и
Т(е) - А ¿17 „ Г)(е) =
а — безразмерная скорость и ; б — коэффициент давления с
1 — первый вариант функции Б0) (х): = 2,53 и = 0,291; 2 — второй вариант функции Б0' (х): = 4,617 и = 0,291;
0
2
4
6
х
0
2
4
6
х
а
На рис. 6 графики зависимостей безразмерной скорости и = и( х) и коэффициента давления с = с (х) от координаты х построены для первого (линии 1) и третьего (линии 2) вариантов выбора функции Б(0)( х).
Анализ данных, представленных на рис. 2-6, позволяет установить взаимосвязь полученных после завершения итерационного процесса конфигураций эквивалентной поверхности с распределениями скорости и давления на обтекаемом теле, а также выявить некоторые особенности, присущие указанным конфигурациям и распределениям. Важным обстоятельством является малое
отличие конечных значений минимального диаметра следа и)' для всех трех рассмотренных вариантов выбора зависимости Б(0) (х) диаметра поперечного сечения участка Е21 эквивалентной поверхности от продольной координаты х . В то же время конечное значение длины следа и окончательная конфигурация хвостовой
части поверхности эквивалентного тела существенным образом зависят от выбора варианта указанной зависимости. Следует отметить,
что если зависимость Б( 0) (х) выбрана и зафиксирована, то дополнительно проведенные расчеты показывают незначительность влияния начальных значений длины и минимального диаметра Д0 следа
на их конечные значения и , а также на окончательную форму эквивалентной поверхности.
0,8
0,4
0,0
Г 2
1
ср 0,8
0,4
0,0
-0,4
1
л
у 2
0 2 4 6 х
4 6 х
а б
Рис. 6. Зависимости параметров потока газа в продольном сечении поверхности эквивалентного тела от координаты х при обтекании тела
с Ь =1 и =2 при сс = 0 и Яе =5-10 :
1 — первый вариант функции Б(0) (х): Де - 2,53 и - 0,291;
I ■ № = '
У«) = I
2 — третий вариант функции Б(х): Д«' - 3,197 и Бу - 0,291; а — безразмерная скорость й; б — коэффициент давления с
Зависимость длины следа и окончательной конфигурации хвостовой части эквивалентной поверхности от выбора варианта функции Б0' (х) приводит к различиям зависимостей безразмерной скорости й (рис. 5а) и коэффициента давления с (рис. 5б) от продольной координаты х . Данные различия сколько-нибудь заметны
0
2
лишь в областях, расположенных вниз по потоку от донного среза, и становятся существенными при х > 4,2 в случае сравнения вариантов 1 и 2 (рис. 5), а также при х > 4,8 в случае сравнения вариантов 1 и 3 (рис. 6).
Таким образом, на поверхности обтекаемого тела (при х < 3) отклонения значений безразмерной скорости и(х) и коэффициента
давления с (х) оказываются несущественными для всех трех
рассмотренных вариантов выбора функции Б^ (х) . Данный факт,
свидетельствует о том, что распределение модуля безразмерной скорости и и коэффициента давления с по поверхности обтекаемого
тела практически не зависит от формы хвостовой части эквивалентной поверхности, которая определяется выбором варианта функции
Б(0) (х). Это дает определенную гибкость в использовании
рассматриваемой методики, так как предоставляет возможность
осуществлять выбор функции Б(0) (х) исходя из специфики решаемой
прикладной задачи. Так, если требуется только найти распределение
модуля безразмерной скорости и и коэффициента давления с по
поверхности обтекаемого тела, то с точки зрения экономии вычислительных затрат рациональным представляется выбор первого или
третьего вариантов функции Б(0) (х). Если же будет решаться еще и
дополнительная задача, например, связанная с определением некоторых характеристик потока в спутном следе, то целесообразным может оказаться выбор второго варианта задания функции Б(0) (х).
Заключение. При использовании методики математического моделирования дозвукового отрывного обтекания удлиненных тел с донным срезом организация итерационного процесса с применением сглаживающих кубических сплайнов позволяет сравнивать распределения скорости и давления по поверхности эквивалентного тела, получаемые для различных вариантов выбора зависимости диаметра поперечного сечения хвостовой части эквивалентной поверхности от продольной координаты. Анализ результатов проведенного вычислительного эксперимента показал малое различие конечных значений минимального диаметра следа при заметном отличии конечных величин длины следа и окончательной конфигурации хвостовой части поверхности эквивалентного тела для рассмотренных вариантов выбора указанной зависимости. Существенно, что данные отличия практически не влияют на распределения скорости и давления по поверхности обтекаемого тела. Это дает определенную гибкость в использовании рассматриваемой
методики, так как предоставляет возможность осуществлять выбор зависимости диаметра поперечного сечения хвостовой части эквивалентной поверхности от продольной координаты исходя из специфики решаемой прикладной задачи.
ЛИТЕРАТУРА
[1] Флетчер К. Вычислительные методы в динамике жидкостей. Т. 2. Методы расчета различных течений. Москва, Мир, 1991, 552 с.
[2] Патанкар С.В. Численные методы решения задач теплообмена и динамики. Москва, Энергоатомиздат, 1984, 152 с.
[3] Гогиш Л.В., Степанов Г.Ю. Отрывные и кавитационные течения. Основные свойства и расчет модели. Москва, Наука, 1990, 384 с.
[4] Белоцерковский С.М., Ништ М.И., Котовский В.Н., Федоров Р.М. Трехмерное отрывное обтекание тел произвольной формы. Москва, ЦАГИ, 2000, 265 с.
[5] Saffman P.G. Vortex Dynamics. Cambridge, Cambridge University Press, 1992, 311 p.
[6] Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Cambridge, Cambridge University Press, 2005, 592 p.
[7] Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. Москва, Ин-т механики МГУ им. М.В. Ломоносова, 2006, 184 с.
[8] Головкин М.А., Головкин В.А., Калявкин В.М. Вопросы вихревой гидромеханики. Москва, Физматлит, 2009, 264 с.
[9] Дергачев С.А., Щеглов Г.А. Моделирование обтекания тел методом вихревых элементов с использованием замкнутых вихревых петель. Научный вестникМГТУГА, 2016, № 223, с. 19-27.
[10] Коцур О.С., Щеглов Г.А. Реализация метода обмена интенсивностями вортонов-отрезков для учета вязкости в методе вихревых элементов. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2018, № 3, с. 48-67.
[11] Kuzmina K.S., Marchevskii I.K., Moreva V.S. Vortex sheet intensity computation in incompressible flow simulation around an airfoil by using vortex methods. Mathematical Models and Computer Simulations, 2018, vol. 10, iss. 3, pp. 276-287.
[12] Тимофеев В.Н. Математическое моделирование отрывного дозвукового обтекания осесимметричных тел с учетом донного давления. Инженерный журнал: наука и инновации, 2014, вып. 10. URL: http://engjournal.ru/catalog/mathmodel/aero/1246.html (дата обращения 17.10.2017).
[13] Тимофеев В. Н. Построение полубесконечного эквивалентного тела при математическом моделировании дозвукового отрывного осесимметричного обтекания. Математическое моделирование и численные методы, 2016, № 4, с. 67-83.
[14] Тимофеев В. Н. Моделирование дозвукового отрывного обтекания тел методом дискретных вихрей на основе концепции эквивалентной поверхности с кубическими сплайнами. Математическое моделирование и численные методы, 2020, № 4, с. 27-43.
[15] Timofeev V.N. Mathematical simulation of the subsonic flow around the lengthening bodies with the flow separation in the region of ground shear with the use of an equivalent body. Journal of Physics: Conference Series, 2018, vol. 1141, art. no. 012095.
[16] Роджерс Д.Ф., Адамс Дж.А. Математические основы машинной графики. Москва, Наука, 2001, 576 с.
[17] Димитриенко Ю.И. Тензорное исчисление. Москва, Высшая школа, 2001, 576 с.
[18] Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент (в математической физике, аэродинамике, теории упругости и дифракции волн). Москва, ТОО «Янус», 1995, 520 с.
[19] Тимофеев В.Н. Моделирование дозвукового отрывного обтекания тел с донным срезом по схеме течения с эквивалентной полубесконечной поверхностью при малых углах атаки. Математическое моделирование и численные методы, 2019, № 4, с. 31-49.
[20] Тимофеев В.Н. Особенности вихревой схемы при моделировании дозвукового отрывного обтекания с полубесконечным эквивалентным телом. Математическое моделирование и численные методы, 2017, № 4, с. 73-91.
[21] Аубакиров Т.О., Белоцерковский С.М., Желанников А.И., Ништ М.И. Нелинейная теория крыла и ее приложения. Алматы, Гылым, 1997, 448 с.
Статья поступила в редакцию 12.09.2021
Ссылку на эту статью просим оформлять следующим образом: Тимофеев В.Н. Методика организации итерационного процесса при моделировании дозвукового отрывного обтекания удлиненных тел с применением эквивалентной поверхности и кубических сплайнов. Математическое моделирование и численные методы, 2021, № 4, с. 80-102.
Тимофеев Валерий Николаевич — канд. техн. наук, доцент кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор более 40 научных и методических работ. Область научных интересов: математическое моделирование, численные методы, механика жидкости и газа, аэродинамика, численные методы, методы оптимизации. e-mail: v_n_1951@mail.ru
Method for organizing an iterative process in modeling subsonic separated flow around elongated bodies using an equivalent surface and cubic splines
© V.N. Timofeev
Bauman Moscow State Technical University, Moscow, 105005, Russia
In the aspect of improving the methodology of mathematical modeling of subsonic detachable flow around elongated bodies with partial implementation of the concept of viscosity-inviolable interaction, the issues of organizing the iterative process of constructing the surface of an equivalent body are considered. A flow diagram with a semi-infinite equivalent surface was used. Numerical modeling was carried out according to the algorithms of the technique using the method of discrete vortices. For greater versatility, when constructing a calculated grid on the surface of an equivalent body, approximation with smoothing cubic splines was used. Data on the influence of the shape of the tail section of the equivalent surface on the distribution of speed and pressure during axisymmetric flow around bodies with a bottom section are presented.
Keywords: mathematical modeling, subsonic separation flow, the concept of visco-inviscid interaction, the equivalent of the surface, smoothing cubic spline
REFERENCES
[1] Fletcher K. Vychislitel'nye metody v dinamike zhidkostej. T. 2. Metody rascheta razlichnyh techenij [Computational methods in fluid dynamics. Vol. 2. Methods for calculating different flows]. Moscow, Mir Publ., 1991, 552 p.
[2] Patankar S.V. Chislennye metody resheniya zadach teploobmena i dinamiki. [Numerical methods for solving problems of heat transfer and dynamics]. Moscow, Energoatomizdat Publ., 1984, 152 p.
[3] Gogish L.V., Stepanov G.Yu. Otryvnye i kavitatsionnye techeniya. Osnovnye svoiystva i raschet modeli [Detached and cavitational flows. Basic properties and model calculation]. Moscow, Nauka Publ., 1990, 384 p.
[4] Belotserkovskiy S.M., Nisht M.I., Kotovskiy V.N., Fedorov R.M. Trekhmernoe otryvnoe obtekanie tel proizvolnoy formy [Three-dimensional detached flow of the bodies of arbitrary form]. Moscow, TsAGI (Central Aerohydrodynamic Institute) Publ., 2000, 265 p.
[5] Saffman P.G. Vortex Dynamics. Cambridge, Cambridge University Press, 1992, 311 p.
[6] Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Cambridge, Cambridge University Press, 2005, 592 p.
[7] Andronov P.R., Guvernyuk S.V., Dynnikova G.Ya. Vikhrevye metody rascheta nestatsionarnykh gidrodinamicheskikh nagruzok [Vortex methods of calculation of nonstationary hydrodynamic loads]. Moscow, Institute of Mechanics Lomon-osov MSU Publ., 2006, 184 p.
[8] Golovkin M.A., Golovkin V.A., Kalyavkin V.M. Voprosy vihrevoj gidromehaniki [Questions of vortex hydromechanics]. Moscow, Fizmatlit Publ., 2009, 264 p.
[9] Dergachev S.A., Scheglov G.A. Vortex element method simulation of flow around bodies using closed vortex loops. Civil Aviation High Technologies, 2016, no. 223, pp. 19-27.
[10] Kocur O.S., Shcheglov G.A. Implementation of the particle strength exchange method for fragmentons to account for viscosity in vortex element method. Herald of the Bauman Moscow State Technical University. Series Natural Sciences, 2018, no. 3, pp. 48-67.
[11] Kuzmina K.S., Marchevskii I.K., Moreva V.S. Vortex sheet intensity computation in incompressible flow simulation around an airfoil by using vortex methods. Mathematical Models and Computer Simulations, 2018, vol. 10, iss. 3, pp. 276-287.
[12] Timofeev V.N. Mathematical modeling of separated subsonic flowaround axially symmetrical bodies with base pressure. Engineering Journal: Science and Innovation, 2014, no. 10. Available at: http://engjournal.ru/ catalog/math-model/aero/1246.html (accessed October 17, 2017).
[13] Timofeev V.N. Construction of a semi-infinite equivalent body in mathematical modeling of subsonic separated axisymmetric flow. Mathematical Modeling and Computational Methods, 2016, no. 4, pp. 67-83.
[14] Timofeev V.N. Modeling of subsonic separation flow around bodies by the method of discrete vortices based on the concept of an equivalent surface with cubic splines. Mathematical Modeling and Computational Methods, 2020, no. 4, pp. 27-43.
[15] Timofeev V.N. Mathematical simulation of the subsonic flow around the lengthening bodies with the flow separation in the region of ground shear with the use of an equivalent body. Journal of Physics: Conference Series, 2018, vol. 1141, art. no. 012095.
[16] Rjgers D.F., Adams J.A. Mathematical elements for computer graphics. New-York, McGraw-Hill Science, 1989, 512 p.
[17] Dimitrienko Yu.I. Tenzornoe ischislenie [Tensor calculus]. Moscow, Vysshaya shkola Publ., 2001, 576 p.
[18] Lifanov I.K. Metod singulyarnykh integralnykh uravneniy i chislennyy eksperi-ment (v matematicheskoi fizike, aerodinamike, teorii uprugosti i difraktsii voln) [The method of singular integral equations and numerical experiment (in mathematical physics, aerodynamics, theory of elasticity and diffraction of waves)]. Moscow, LLP Yanus Publ., 1995, 520 p.
[19] Timofeev V.N. Simulation of the subsonic detachable body with a bottom cut on the current pattern with an equivalent half-infinite surface at small angles of attack. Mathematical Modeling and Computational Methods, 2019, no. 4, pp. 31-49.
[20] Timofeev V.N. Special features of vortex diagram in simulation of subsonic detached flow around the semi-infinite equivalent body. Mathematical Modeling and Computational Methods, 2017, no. 4, pp. 73-91.
[21] Aubakirov T.O., Belotserkovskiy S.M., Zhelannikov A.I., Nisht M.I. Nelineinaya teoriya kryla i ee prilozheniya [Nonlinear wing theory and its applications]. Almaty, Gylym Publ., 1997, 448 p.
Timofeev V.N., Cand. Sc. (Eng.), Assoc. Professor of Computational Mathematics and Mathematical Physics Department, Bauman Moscow State Technical University. Author of over 40 research and methodological papers. Research interests include mathematical modeling, computational methods, fluid mechanics, aerodynamics, optimization methods. e-mail: v_n_1951@mail.ru
