Расчет пространственного обтекания тела вращения в сжимаемом потоке газа методом гидродинамических особенностей

Расчет пространственного обтекания тела вращения в сжимаемом потоке газа методом гидродинамических особенностей

П.И. Гайда, к.т.н.;

С.Д. Косторной, проф., д.т.н.

(Сумский государственный университет)

Рассмотрена нестационарная нелинейная задача о произвольном движении тела вращения в сжимаемом газе. Возмущенное телом движение газа считается потенциальным. Решение задачи ищется в виде суммы запаздывающих потенциалов диполей, распределенных по поверхности тела, и сходящей с его задней кромки вихревой пелене, и источников, распределенных вне тела и следа. С помощью источников учитываются нелинейные члены в уравнении для потенциала возмущения скорости. Численный метод строится на основе перехода от непрерывных распределений особенностей по поверхности тела к панельным и к дискретным во времени. Вихревой след на теле и в следе заменяется системой рамок диполей с равномерной интенсивностью по поверхности рамки. Интенсивности рамок диполей определяются последовательно, шаг за шагом по времени. В процессе расчета выстраивается вихревая пелена и определяются интенсивности источников.

Рассмотрим неустановившееся обтекание тела идеальным сжимаемым газом, в котором отсутствуют процессы передачи тепла, обусловленные свойством вязкости (теплопроводность, диффузия), и газ не излучает энергию, которое в аэродинамике называют изоэнтропическим (адиабатическим). Предполагается, что массовые силы отсутствуют и нет скачков уплотнения, а завихренность потока сосредоточена внутри тонких слоев, сходящих с кромок. Эти слои моделируем тонкими поверхностями, на которых терпит разрыв касательная к поверхности - составляющая скорости газа.

Состояние газа будет описываться уравнениями неразрывности, движения и постоянства энтропии:

(1.1)

Здесь - скорость движения газа; - плотность; - давление; - показатель адиабаты.

Считая течение вне поверхностей свободных вихрей потенциальным, т.е. , введем потенциал скорости Ф. Тогда

, где- оператор градиента, (1.2)

а преобразованное уравнение неразрывности будет представлять собой уравнение для потенциала скорости Ф [1]:

(1.3)

где производные

В векторной форме уравнение (1.3) может быть записано следующим образом [2]:

. (1.4)

Здесь а2 - квадрат скорости звука, равный

, (1.5)

показатель адиабаты.

Следствием уравнений движения является интеграл Коши - Лагранжа [1]:

, (1.6)

из которого получается выражение для давления

, (1.7)

в котором индексом «» обозначены параметры невозмущенного течения.

Отметим, что из уравнения (1.1) также следует, что сошедшие в поток свободные вихри должны двигаться вместе с частицами жидкости.

Граничные условия, необходимые для решения задачи обтекания тела, будут иметь вид:

1) условие непротекания на теле

,(1.8)

где нормаль к поверхности тела;

- вектор скорости движения ее точек;

2) на вихревых поверхностях свободных вихрей , сходящих с кромок (рис. 1), должно выполняться условие непрерывности нормальной составляющей скорости и давления

.(1.9)

Последнее условие может быть записано в виде

(1.10)

где [Ф] - скачок потенциала на поверхности,

скорость на вихревой поверхности.

Рисунок 1 - Расчетная схема обтекания тела вращения

Из формулы (1.10) следует, как уже было отмечено, что точка свободной вихревой пелены должна перемещаться вместе с жидкими частицами. Необходимо также потребовать, чтобы

,(1.11)

а также задать начальные условия. Для газа, покоящегося при t=0, оно будет иметь вид

(1.12)

Для определения положения вихревой пелены используется уравнение движения жидкой частицы

.(1.13)

При решении задачи обтекания тела уравнение (1. 4), записанное в виде

(2.1)

нужно решать шаг за шагом по времени, считая на каждом последующем шаге правую часть уравнения (2.1) известной из решения на предыдущем шаге. При начальных условиях (1.12) в момент t=0 ее можно считать равной нулю. Решение волнового уравнения (2.1) можно записать в виде [3]:

(2.2)

где [Ф]=Ф+ - Ф- - скачок потенциала на теле S и следе ,

(t) - дельта-функция Дирака.

Скорость возмущенного потока определяется дифференцированием (2.2):

.(2.3)

Преобразуем первый интеграл, соответствующий поверхностному распределению диполей. Дифференцируя его, получим

.(2.4)

С помощью соотношения

,(2.5)

справедливого, если вектор не зависит от , преобразуем формулу (2.4) к виду

.(2.6)

Можно убедиться, что

.(2.7)

Оператор касателен в каждой точке к поверхностям S и . Учитывая, что

,(2.8)

и используя аналог формулы Стокса, из формулы (2.6) получим

(2.9)

,

где контур L обходится против часовой стрелки, если смотреть со стороны нормали.

Вводя вектор интенсивности вихревого слоя

,(2.10)

Получим

(2.11)

.

Второй интеграл в этой формуле можно трактовать как вихрь, сосредоточенный на границе области, третий соответствует вихрям, распределенным на поверхностях S и .

Процесс изменения интенсивностей рамок диполей на теле и в следе дискретизируется во времени. Непрерывное во времени изменение интенсивностей диполей заменяется кусочно - постоянным. Будем считать, что скачкообразное изменение интенсивностей диполей на теле происходит непосредственно за моментами времени s (S=0,1,…,-1), где 0=0. Граничное условие (1.8) на теле будем удовлетворять в моменты r (r=1,2,…,) в центре каждой из рамок. Считаем, что на интервале времени |0,t1| вихревая пелена отсутствует, а также отсутствуют, в силу начальных условий (1.11), источники в газе. Вычисляя скорости от всех рамок диполей в расчетных точках на поверхности тела и удовлетворяя граничным условиям, получаем систему линейных уравнений, из решения которой определяются интенсивности диполей на теле в первый момент времени. Эти интенсивности используются затем для вычисления возмущенных скоростей на кромке тела, что необходимо для определения движения рамок диполей в следе.

Процесс формирования следа описывается движением и деформацией рамок диполей следа. Считается, что примыкающая к кромке тела сторона этой рамки движется со скоростью кромки, сторона рамки, сходящая в след, - со скоростью газа. Интенсивности диполей этой рамки принимаются равными интенсивностям диполей на кромке тела на предшествующем интервале времени. Скорости от рамок диполей на теле и в следе используются затем для вычисления интенсивностей источников в возмущенном телом объеме газа.

В последующем для вычисления интенсивностей диполей на крыле и за ним используется след, найденный на предыдущих временных шагах. Так, к моменту 2 будут известны интенсивности рамок диполей в следе, их форма и положение, а также интенсивности объемных источников в возмущенной телом области. Считая, что на теле интенсивности диполей изменялись в момент 1, можем найти интенсивности диполей на теле в момент 2, выполнив граничное условие (1.8) на теле. Затем вычисляются скорости движения рамок в следе, и численным интегрированием уравнения (1.13) определяется новое положение рамок диполей. При этом интенсивности рамок, образующихся у кромки тела, принимаются равными интенсивностям диполей на кромке тела, интенсивности уже образовавшихся рамок известны. После этого процесс повторяется для последующих расчетных моментов.

Скорости, индуцируемые плоской четырехугольной рамкой диполей в сжимаемом газе, получены С.Д. Шипиловым в [5], а одинаковые с ними - от замкнутого суммарного нестационарного вихря в дозвуковом сжимаемом потоке при скачкообразном изменении напряженности - В.И. Бушуевым в работе [6].

Для рамки диполей выражение для скоростей имеет вид

.(3.1)

Здесь циркуляция по контуру, проходящему через рамку в момент ; вектор, характеризующий положение i-й вершины рамки в момент ; t - время, проходящее после образования рамки;

Вид функции

приведен в работе [5]. Как там показано, при дозвуковых скоростях движения рамки функция WD отлична от нуля только при (0 - время испускания сигнала плоскостью рамки, пришедшего в момент t в точку ). Она не зависит от t при (1- максимальное значение , при котором в точке в момент t «слышимы» все вершины рамки). Поэтому интеграл в (3.1) определяется численно лишь в достаточно небольшом диапазоне .

Представив интеграл в формуле (3.1) в виде

,(3.2)

и сделав замену переменных

,

получим, считая на интервале постоянной,

.

Для вычисления функции на малом интервале используется линейная аппроксимация закона движения панели

.(3.3)

Такая кусочно - линейная аппроксимация закона движения угловых точек панели позволяет с достаточной точностью учесть запаздывание прихода сигнала в точку . В результате выражение для скорости может быть представлено в виде

i=1,…,4. (3.4)

Для определения интенсивностей диполей на теле введем следующую систему обозначений для координат характерных точек на теле и в следе. Пусть n - число разбиений тела на панели вдоль характерного размера (длина тела вращения); N - число разбиений тела вдоль окружности сечения тела. Будем отмечать индексами i, k левый верхний угол на нем, индексом S - момент времени. Обозначим через вектор координат этой точки в момент S в неподвижной системе координат, через - циркуляцию панели с номером (i, k) в момент S.

Нумерацию панелей в следе будем вести в соответствии с номером того расчетного момента, в который эта панель начала образовываться, и в соответствии с нумерацией панелей вдоль той кромки, с которой панель сошла. Нумерацию панелей в следе будем отмечать индексом . Значение индекса =1 означает, например, что панель начала сходить в момент 1. Обозначим через вектор координат левой верхней точки панели с номером (, k) в следе за задней кромкой, через - вектор координат левой верхней угловой точки в момент . Аналогично обозначим и интенсивности панелей в следе и . Интенсивности диполей этих панелей равны, в силу сохранения циркуляции, интенсивностям соответствующих панелей диполей тела, примыкающих к кромке, в момент схода (т. е. в момент Н), поэтому можем записать

(3.5)

При составлении системы уравнений для определения интенсивностей диполей на теле в момент r к этому времени должны быть известными интенсивности, положения диполей в следе и интенсивности источников в возмущенном телом объеме. Выполнив условие непротекания на теле в расчетных точках, лежащих в центрах панелей тела, и вычислив в них скорости от всех панелей с учетом запаздывания возмущений в расчетных точках, получим

(3.6)

.

Здесь вектор скорости, индуцированной панелью диполей единичной интенсивности за интервал времени |s-1-s| в расчетной точке с номером j.

Вектор можно представить в виде [5]:

+,(3.7)

где если не принадлежит интервалу [];

если принадлежит интервалу [];(3.8)

если не принадлежит интервалу [];

если принадлежит интервалу [].(3.9)

Входящие в формулу (3.7) и равны соответственно:

,(3.10)

,

а и следует вычислять по формуле

.(3.11)

В формуле (3.11) использованы следующие обозначения:

,

,(3.12)

.

, , .

Изложенный алгоритм расчета нестационарного обтекания тела вращения 3Д потоком сжимаемой жидкости реализован на ЭВМ.

Рисунок 2 - Векторное поле скоростей в окрестности тела вращения при угле атаки 450 за 10 моментов времени

Для тела (рис. 1) диаметром 70 мм, длиной l=300 мм решена прямая задача для случая V?=100м/с и различных углов атаки. На рис. 2 показано векторное поле скоростей в окрестности тела с учетом вихревого следа за 10 моментов времени. Зависимость аэродинамических коэффициентов Сх, Су, См в начальный момент времени от угла атаки представлена на рис. 3.

Рисунок 3 - Зависимость аэродинамических показателей тела вращения от угла атаки в начальный момент времени

Выполненные расчеты и численные значения коэффициентов подтвердили возможность численного решения трехмерной задачи нестационарного обтекания тел сжимаемым потоком газа методом гидродинамических особенностей.

СПИСОК ЛИТЕРАТУРЫ

Белоцерковский С.М., Скрипач Б.К., Табачников В.Г. Крыло в нестационарном потоке газа. - М.: Наука, 1971.

Бушуев В.И. К расчету пространственного поля скоростей нестационарного вихря неизменной формы при дозвуковых скоростях // Научно-методические материалы по аэродинамике летательных аппаратов. ВВИА им. Н.Е. Жуковского. - 1987. -С.192-202.

Владимиров В.С. Уравнения математической физики. - М.: Наука, 1967.

Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. - М.: Физматгиз, 1963. - Ч.1.

Сирс У.Р. Общая теория аэродинамики больших скоростей. - М.: Изд-во МО СССР, 1962.

Шипилов С.Д. Расчет возмущений от панели диполей, произвольно движущейся в сжимаемом газе // Труды ВВИА им. проф. Н.Е. Жуковского. - 1986. - Вып. 1313. - С.488-500.