Усовершенствование метода расчета трансзвукового обтекания тел вращения Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м VI 1975

№ 6

УДК 533.6.011.3

УСОВЕРШЕНСТВОВАНИЕ МЕТОДА РАСЧЕТА ТРАНСЗВУКОВОГО ОБТЕКАНИЯ ТЕЛ ВРАЩЕНИЯ

Р. К. Тагиров

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

1. Эффективный метод расчета обтекания тел вращения типа мотогондол или фюзеляжей дозвуковым или трансзвуковым потоком был развит в работах [1, 2]. В методе используется разностная схема С. К. Годунова [3, 4] и решение достигается в процессе установления по времени. Метод позволяет проводить расчеты при наличии втекающего в двигатель и истекающего из двигателя потоков и определять распределение параметров в поле течения и величины сопротивления. Расчет с помощью метода установления проводится для идеального (невязкого и нетеплопроводного) газа. Различные реальные эффекты (пограничный слой, донный уступ, эжектирующее действие струи и т. д.) учитываются с помощью дополнительных приближенных методов.

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

1) Уточнение граничного условия на внешней границе. В работах [1, 2] в качестве граничного условия на внешней границе области течения бралось или равенство нулю нормальной составляющей скорости (VB — 0), что имитировало непроницаемую стенку, или условие постоянства статического давления {рв = const), что имитировало идеально проницаемую стенку трубы. Необходимый минимальный радиус внешней границы, при котором отсутствует

влияние границы на распределение параметров вдоль тела, устанавливался на основе параметрических расчетов для разных радиусов. Следует отметить, что граничные условия рв — const, и тем более VB = 0, не соответствуют распределению параметров вдоль внешней границы, которое возникает при обтекании тела безграничным потоком. В то же время ясно, что начиная с некоторого радиуса внешней границы возмущения, приходящие от тела на эту границу, будут малы и, следовательно, для определения параметров на этой внешней границе можно использовать соотношения линейной теории, как это было предложено в работе [5|.

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

(1 _М2 ) д2 9 4- _ (х+ ^ М” дф д*<?

00 дх2 ду2 U со дх дх2

(х + 1) М2, df .. , „

При --------------- —- — Cl это уравнение можно считать линей-

ным. Следовательно, если возмущения, приходящие на внешнюю границу от тела, таковы, что выполняется условие &<1, то можно воспользоваться решением, полученным в рамках линейной теории. Ограничимся рассмотрением случая дозвукового течения. Составляющие скорости возмущения определяются следующим образом:

а, = ^ г __ /(£) {х - S)dg_

дх J [(х-_|)3 + (1_м^)Я3/2 ’

О

V' = _*L = (1 _ М2 )у Г-----------------------;

ду 00 J [(jf_6)i + (l_ М^)у2]3'2

о

здесь / — длина рассматриваемой области течения; у=ув = const — радиус внешней границы. Функции /($) пока неизвестны.

Полученные соотношения для и! и v' необходимо связать с соотношениями, используемыми в методе установления. Для этого используется следующий итерационный алгоритм. Вначале принимается /(£)еееО, что эквивалентно принятию граничного условия рв = const, и проводится расчет обтекания методом установления. Через z временных шагов в поле течения вырабатываются промежуточные значения параметров, с помощью которых можно определить, применяя линейную интерполяцию по у, распределение вертикальной составляющей скорости vB — v' вдоль внешней границы области течения. Выделим на этой границе я точек, в каждой из них известны значения v. (j меняется от 0 до га). Применяя формулу Симпсона для аппроксимации интеграла в выражении для v', можно получить систему уравнений:

____________/о____________ J_ ___________4Л___________ , _ '

[(*; - ?о)2 + (1 - Ml) у\ ]3/2 ‘ [(*,• - Si)2 + (1 - л£) У2В ]3'2 ' ' '

_j____________fn____________ Зя У,

’’’ [(^--У2+(1-М / (1-м 1)ув '

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

выделенной точки значения и'.:

'' Зп f [(-*>• — So)2 Ч- (1 -м^)^]3

/п (Xj Sn)

и соответствующие значения давления Р) = Р

М2 / а-2 4- v - 2

здесь Moo, t/oo, Poo и х — число М, скорость, давление и показатель адиабаты невозмущенного набегающего потока.

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

Следует отметить, что при уменьшении ув возмущения, приходящие на внешнюю границу от тела, возрастают, поэтому минимально допустимое значение ув ограничивается написанным выше условием &<С1- Выполнение этого условия, гарантирующего применимость результатов линейной теории, всегда проверялось при проведении расчетов. Для установления сходимости решения и уменьшения времени счета при использовании нового граничного условия на верхней границе (pB = var) были проведены методические расчеты на ЭЦВМ М-222. Исследовалось обтекание кормовой части тела вращения, образующая которой представляла собой дугу окружности радиуса Л? = 9,68 (за единицу длины принят максимальный радиус тела). Слева и справа к кормовой части примыкали цилиндры. Расчеты методом установления велись до достижения заданной точности. Все примеры расчета получены при п= 10 и х= 1,4. Если не делается оговорок, то в качестве начальных параметров для каждого примера расчета используются данные, полученные из расчета одномерного изэнтропического потока.

Был проведен расчет обтекания при использовании старого граничного условия (рв = const) и разностной сетки (30X10). показанной на фиг. 1 (Моо — 0,9, ув = 5). Полученные распределения коэффициента давления вдоль тела (ср„) и вдоль верхней границы расчетной области (срв) показаны на фиг. 2 пунктирной линией. При этом получен коэффициент волнового сопротивления схв = =—0,0005^0 и истрачено машинного времени —7 ч. Результаты расчета обтекания этого же тела при использовании уточненного граничного условия (рв — маг) показаны на фиг. 2 штриховой линией для j/B = 3 и сплошной линией для ув = 2. В связи с более точным граничным условием были взяты меньшие значения ув и в два раза меньшее количество расчетных ячеек в направлении оси_у. Соответствующие разностные сетки (30X5) показаны на фиг. 1.

----1-----1--1----------1----1___I____I I_____I ,

О 1 2 3 4 5 6 7 8 9 10 X

Фиг. 2

Из этих расчетов найдены значения схв =—0,003 ^0, к — 0,22 для ув = 3 и схв= — 0,01, /2 = 0,38 для ув = 2. Причем для расчета первого варианта (^/в = 3) было истрачено 2 ч 40 мин, для второ-

го — 3 ч машинного времени.

Из фиг. 2 видно, что распределения срп полученные при старом и новом граничных условиях, практически совпадают, хотя в последнем случае было затрачено в 2,5 раза меньше машинного времени, чем в первом случае. Результаты для ув = 2 получились

менее точными, так как из-за возрастания приходящих от тела возмущений (максимальное значение & = 0,38) использование результатов линейной теории становится менее обоснованным.

Следует отметить некоторые методические вопросы получения решения. Вначале была сделана попытка получения решения при Моо = 0,9, когда поправка граничного давления рв проводилась через 2 = 50 временных шагов. Получить установление по времени в этом случае не удалось. Однако после принятия z =100 решение было получено. Используя полученные результаты в качестве начальных данных, были проведены расчеты для 2 = 50 и 10. В том и другом случае было достигнуто установлен-ие по времени, результаты совпали с данными, полученными при z == 100.

Для получения представления о сходимости решения при изменении числа Mo, были проведены расчеты обтекания рассматриваемого тела при Моо = 0,94, ув = 5 и 10. При использовании нового граничного условия и сетки (30X5) для обоих вариантов решение было получено при z = 200. При этих расчетах найдены значения: схв = 0,0132, й = 0,22 при ув = 5 и с*в = 0,0138, 6 = 0,06 при ув — 10. На каждый из вариантов было истрачено приблизительно 3 ч машинного времени.

2) Использование переменного шага по времени по ячейкам. При проведении расчетов методом установления шаг по времени должен быть достаточно мал, чтобы волны, возникающие при распаде произвольного разрыва, не достигали соответствующих противоположных граней ячейки [3]. В работах [1, 2.] для определения шага по времени т использовался следующий алгоритм. После определения в некоторый момент времени t параметров потока во всем поле течения для каждой ячейки определялись: ъх—время, потребное для достижения волной противоположной грани ячейки в направлении оси х; — время, потребное для достижения волной противоположной грани ячейки в направлении, нормальном к

продольной грани ячейки; т = —----------суммарный шаг по време-

ZJC "Ь Zy

ни [3]. Из всех значений т, полученных в ячейках, определялось минимальное значение xmin, которое и принималось в качестве последующего шага по времени.

Однако, если не определять тга)П, а расчеты для каждой ячейки вести со своим индивидуальным шагом то можно получить решение за меньшее время. С целью проверки этого предположения были проведены расчеты обтекания кормовой части, рассмотренной выше, при Моо = 0,9. При использовании сетки (30 X Ю), показанной на фиг. 1, старого граничного условия рв = const и индивидуального шага по времени т для каждой ячейки на решение задачи было истрачено 5 ч машинного времени. Причем найденные результаты, в том числе схв = — 0,0002 :=0, очень точно совпали с результатами для исходного варианта (пунктирная линия на фиг. 2), полученными после 7 ч счета на ЭЦВМ. Таким образом, видно, что использование в каждой ячейке своего шага по времени позволяет приблизительно в 1,4 раза сократить время счета.

Интересно определить суммарное сокращение времени счета при применении и нового граничного условия, и переменного шага по времени по ячейкам. Такой расчет был проведен для сетди <30 X 5), показанной на фиг. 1 для ув = 3. Результаты расчета очень

точно совпали с результатами, полученными при использовании лишь нового граничного условия (штриховая линия на фиг. 2), для получения решения было израсходовано 2 ч 20 мин машинного времени.

3) Уточнение начальных данных. Для всех рассмотренных выше примеров в качестве начальных данных использовались параметры, найденные из расчета одномерного изэнтропического течения. Если эти начальные данные определять точнее, то можно дополнительно сократить время счета. В связи с этим был разработан приближенный метод расчета, основанный на методе работ [6, 7]. В последней из указанных работ для скорости возмущения на поверхности тела получено соотношение

где

$ХХ (•*)

1п X -4-

и..

I

тг

Х=1-М»,-

М.2

* +1

00 и„

8л:

I

[$(?)

2 (лс — 6) +

\{х-

.Ц)2 + ±8(Х)

К

, аМ

3/2

йх <1Хъ

(5 —площадь поперечного сечения тела).

Для получения решения написанного уравнения в работах [6, 7] дифференцируют его, а затем проводят численное интегрирование полученного дифференциального уравнения. Однако наличие второй и третьей производных от площади поперечного сечения в дифференциальном уравнении создает значительные трудности при расчете обтекания используемых на практике конфигурации. В связи с этим в данной работе рассматривается приближенное решение написанного выше уравнения. Это тем более оправдано, что найденное решение будет использоваться лишь для получения начальных данных для метода установления.

Будем считать, что Моо<1, тогда 1пХ можно разложить в ряд и отбросить члены, содержащие и,' в квадрате и в более высоких степенях. После подстановки в исходное уравнение и проведения обезразмеривания получим следующее решение в явном виде:

где

•ЗД = 2ун

Лук

<14

и„

а $хх (-*■)

ио.

1

2

$(■*)-.У*. 5,(х) = 2У*

2 + Ь8ХХ(*)

{х-Ц) + Зх(х) [(^ — 5)2 + 5 (лг)]3'2

Луя

(1х

■Л,

а = М2 (1 + —М2 + —Ь 00 \ 2 00 3 ” 1 /

(* + 1)М2

Составляющие вектора скорости ин и ■»„, а также коэффициент давления на поверхности тела могут быть определены из следующих соотношений:

мн — 1/со -Т мн

с1х

<1ун

йх

При рассмотрении некоторых тел возникает сложность из-за наличия в выражении для и'н второй производной от площади поперечного сечения, так как на поверхности тела всегда имеются точки стыковки и излома. Для избежания этой сложности в данной

„ № Ун

работе проводилось сглаживание производной —— и вводился по-

йх2

правочный коэффициент при определении Бхх(х) из условия совпадения получаемых ср с результатами метода установления.

Таким образом, для приближенных расчетов использовалось следующее выражение для второй производной:

Для определения параметров в поле течения при известных параметрах на поверхности тела предполагалось, чго распределения и! и V’ описываются экспоненциальным законом:

и' — инеун у.

V'

еу“ у

В этом случае все искомые параметры в каждой ячейке поля течения определяются с помощью следующих соотношений:

« = £/„ + и', Ч) — !)’,

Г . . (и')2 +(„')* 1^

Р = Рсо 1— (х— 1)

Для сравнения были проведены расчеты обтекания рассмотренной ранее кормовой части с помощью приближенного метода и метода установления. Полученные величины давления р в поле течения показаны на фиг. 3. Цифры 1, 2, 3, 4, 5 означают номера

Р

0,96

0,92

ОМ

ом

0,80

0,76

0,72

М„=0,9> ув*3

Г

* “ 5^

З^Л о //! 7/

/ V

Фиг. 3 0 1 2 3 Ч 5 6 7 8 9 х

ячеек, отсчитываемые от поверхности тела в поперечном направлении. На фиг. 3 для сравнения пунктирной линией показано распределение давления, получаемое в случае одномерного течения. Следует отметить, что для расчета методом установления при использовании нового граничного условия, переменного по ячейкам шага по времени и результатов, полученных с помощью приближенного метода, в качестве начальных данных потребовалось 2 ч 10 мин машинного времени, т. е. достигнуто суммарное сокращение времени счета на ЭЦВМ более чем в три раза.

2. При использовании метода установления все поле течения разбивается на четырехугольные ячейки, к каждой из них применяются разностные уравнения [3, 4], позволяющие по известным параметрам в момент времени I определить параметры потока в момент времени t-\-1 (явная разностная схема с первым порядком аппроксимации). Для повышения точности расчетов в местах, где ожидаются относительно большие градиенты параметров, ячейки берутся меньших размеров, а в областях, где ожидаются относительно малые градиенты, ячейки берутся больших размеров, т. е. поле течения разбивается на ячейки разных размеров с переменным шагом. Однако опыт расчетов показывает, что при этом не всегда достигается ожидаемое повышение точности расчета. Например, был проведен расчет обтекания кормовой части, образующая которой описана дугой окружности радиуса 6,3, при Моо =0,7 й ув = 5. Была взята неравномерная разностная сетка (30 X Ю). Найденная в этом случае величина коэффициента волнового сопротивления сЛ.в^0,028 значительно отличается от точного значения схв= = 0. Когда была взята более равномерная сетка, то оказалось, что схв = 0,004 ^0, т. е. получено почти точное значение.

Итак, разностная сетка может влиять на результаты расчета, причем более равномерная сетка позволяет находить более правильные значения схв. Однако при ограниченном количестве ячеек в поле течения (ограниченные память ЭЦВМ и машинное время) брать слишком равномерную сетку нежелательно, так как в этом случае погрешности могут возрасти из-за наличия более крупных ячеек в областях с большими градиентами параметров течения.

Рассматриваемое течение описывается дифференциальными уравнениями в частных производных, которые при решении задачи записываются в виде соответствующих разностных уравнений. Последние должны аппроксимировать исходные уравнения с первым порядком точности (разностная схема С. К. Годунова [3]). Однако при наличии неравномерной сетки возникают погрешности аппроксимации. Для подтверждения рассмотрим представление в разностном виде дифференциального уравнения, причем для простоты ограничимся рассмотрением уравнения неразрывности для одномерного течения

др , д?и _ д ,

дЬ дх

Разобьем поле течения на ячейки, и для некоторой ячейки, ограниченной точками Х-, и лгг_!., запишем соответствующее разностное уравнение :

где х — шаг по времени, кг = х1 — х,_1, /? и и — плотность и скорость на соответствующих границах ячейки.

В соответствии со схемой С. К. Годунова величины /? и £/ определяются из рассмотрения распада произвольного разрыва [3, 4]. Если скорости обеих волн, возникающих при распаде, положительны (сверхзвуковое течение), то как известно [3, 4], параметры Я и и берутся равными значениям р и и в соответствующих левых ячейках, т. е. (Ии)1 = (ри)[ и (#60г-1 = (ри)^_1. Разностное уравнение примет теперь вид

Р;, <+т-Р/,< , (Р«'/- (Р«)г-1 ■ Л

-----------1--------Г----= и .

•с Л;

Разложения в ряд Тейлора функций (р«),_1 и рг, <+- дают:

/ дои \ Л; + Л; ! , ,

(р«)/-, = (Р«),-[-£-1 + о ((А, + А,_о2) ;

' Р«, * + х = Р/, ( + (-^). х +■ О (т2) .

Подстановка в разностное уравнение приводит его к виду

№ + Л/-])2

-ш + °М+ лг-ай;—

= 0.

Видно, что если Лг = А/_ь то написанное выше разностное уравнение аппроксимирует исходное дифференциальное уравнение с первым порядком точности. Однако при такой точности

аппроксимации уже не получается, так как величина — 1^'~1 отличается от единицы. В связи с этим и возникают погрешности при использовании неравномерной сетки в расчетах. Для исключения этих погрешностей необходимо модернизировать разностную схему. Пусть имеются две ячейки, следующие друг за другом и имеющие длины к1 и В середине каждой из этих ячеек из-

вестны параметры р„ и1 и р,-_ь ш—\. Введем основное предположение: при рассмотрении распада разрыва между большей и меньшей ячейками параметры для большой ячейки берутся не в середине ячейки, а на расстоянии, равном половине длины меньшей ячейки, причем сами эти параметры определяются линейной интерполяцией по известным параметрам в серединах двух рассматриваемых ячеек, т. е., если например, Аг_ то разностное уравнение будет

Р/.Н-т —Р 1,1 , (рм)г — (рц),-_1 _п

х Лг '

где черта сверху означает, что параметры ячейки г — 1 берутся на расстоянии /гг/2 от границы, причем

— 9/?-

(р«),_ 1 = (ри)1 + [(ри)/_1 - (ри),] -щ-+ ^ .

Отметим, что аналогичная линейная интерполяция для определения распределения функции внутри ячейки используется в работе [8], посвященной повышению порядка точности в случае равномерной сетки. Теперь можно проверить порядок аппроксимации.

После подстановки полученного выражения для (ри)г-і во второй член уравнения и разложения в ряд Тейлора получаем

' тт^-1 -1£- + о

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

З I—[ "1 -1- І ггп I 11 11 I I I ГI I I 1 I I I I 1 1 1

О_____і____і____і_:______і____і——і.

Фиг. 4

Фиг, 5

С целью проверки предложенного уточнения были проведены специальные расчетные исследования. Было рассмотрено обтекание кормовой части, образующая которой представляла собой дугу окружности радиуса = 9,68 при Моо = 0,9, для двух разностных сеток. В одном случае была взята относительно равномерная разностная сетка, а во втором случае сильно неравномерная (фиг. 4). В обоих случаях было принято ув = 3 и использовалось новое граничное условие (pв ■■=var). Результаты расчета для относительно равномерной сетки при использовании старой разностной схемы и уточненной разностной схемы показаны на фиг. 5 сплошной и штриховой линиями соответственно. Как и следовало ожидать, эти кривые распределения ср вдоль тела близки друг к другу. Соответствующие результаты расчета для сильно неравномерной сетки показаны на фиг. 5 штрих-пунктирной и пунктирной линиями. Видно, что эти распределения существенно отличаются друг от друга. При этом распределение ср (пунктирная кривая), полученное с использованием уточненной разностной схемы, удовлетворительно согласуется с результатами, полученными для относительной равномерной сетки, в то время как использование старой разностной схемы дает искаженный результат.

В заключение автор благодарит А. Н. Крайко за полезные советы.

ЛИТЕРАТУРА

1. Тагиров Р. К. Расчет обтекания кормовых частей тел вращения дозвуковым или трансзвуковым потоком. „Изв. АН СССР, МЖГ“, 1972, № 6.

2. Тагиров Р. К. Трансзвуковое обтекание тела вращения при истечении реактивной струи из его кормовой части. „Изв. АН СССР, МЖГ\ 1974, № 2.

3. Годунове. К., Забродин А. В., Прокопов Г. П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной. „Журн. вычисл. матем. и мат. физ“, ^ I, № б, 1960.

4. И в а н о в М. Я., К р а й к о А. Н. Численное решение прямой задачи о смешанном течении в соплах. „Изв. АН СССР, МЖГ“ 1969, № 5.

5. К р а й к о А. Н., Тагиров Р. К. К околозвуковому обтеканию тел вращения с протоком при наличии истекающей из протока струи. „Изв. АН СССР, МЖГ“, 1974, № 6.

6. Spreiter J. R., Alksne A. Y. Slender-body theory based on approximate solution of the transonic flow equation. NASA Technical Report R-2, 1959.

7. Presz W., Konarski М., Grund E. Prediction of installed nozzle flow fit Ids. AIAA Paper N 70-700, 1970.

8. Колган В. П. Применение принципа минимальных значений производной к построению конечноразностных схем для расчета разрывных решений газовой динамики. „Ученые записки ЦАГИ-, т. III, № 6, 1972.

Рукопись поступила 8jX 1974 г.