СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИЙ МЕТОД ДЛЯ НЕСТАЦИОНАРНЫХ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ

1. Следуя [115], выпишем в обычной тензорной записи дифференциальные законы сохранения массы, импульса и полной энергии для однокомпонентной сплошной среды

которые не зависят от конкретной модели среды и являются незамкнутой системой. Здесь р — плотность, V - вектор относительной скорости с физическими компонентами и v2, Уз» °kj — физические компоненты тензора напряжений; Е - е + V2/2 - полная энергия единицы массы, е = е(р, Т) — удельная (на единицу массы) внутренняя энергия; Q - объемная плотность источников (стоков) тепла; W = | wl9 w2, w3 j = Krgrad T — поток тепла за счет теплопроводности для изотропной среды, — коэффициент теплопроводности, Т — температура; в правых частях уравнений импульса и энергии включают возможные массовые силы, члены, связанные с криволинейностью принятой ортогональной системы координат хх, х2, х3, а в случае, если система координат xlf х2, х3 неинерциальна, также члены, связанные с переносной скоростью движения этой системы относительно некоторой инерциальной системы координат. В частности, в инерциальной декартовой системе координат х, у, z при отсутствии массовых сил /у = fE = 0.

Для замыкания системы (1) —(3), как известно, необходимо привлечь определяющие соотношения, характеризующие конкретную среду. В случае невязкого, нетеплопроводного газа, который рассматривается в данной главе, такими соотношениями являются: предположение о сферичности тензора напряжений

уравнение состояния, связывающее, например, давление р, плотность р и внутреннюю энергию е

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

В случае необходимости используются также калорическое уравнение состояния е = е(р, Т) и другие термодинамические соотношения. В данной главе в качестве одной из искомых функций выбирается удельная энтальпия И = е + р/рив соответствии с этим, уравнение состояния (5) используется в виде

в частности, для совершенного газа

где к = cp/cv = const - показатель адиабаты.

170

Выбирая в качестве искомых величин вектор и = р9 v v2, v3, h , запишем трехмерные нестационарные уравнения течения газа (1) —(7) в виде (1.2.2) из гл. 1

в частности, в случае совершенного газа (8)

Вектор правых частей f = 1 /о>/1*/2»/з> АI» как отмечалось, зависит от выбора системы координат и при отсутствии массовых сил для инерциальных ортогональных криволинейных систем координат с коэффициентами Лямэ Н1. И,. Н* имеет вид

2. При сверхзвуковом пространственном обтекании достаточно гладких затупленных тел конечных размеров или ограниченных в поперечном направлении потоком газа, параметры которого рж, , рж, и т.д. являются известными функциями точки г и времени /, наблюдается картина течения с отошедшей ударной волной, отделяющей область, в которой проявляется влияние тела (’’возмущенная область” или ударный слой) от всего остального потока (’’невозмущенная область”) (рис. 5.1). Эта примыкающая к телу область возмущенного течения может изменять со временем свою форму вследствие движения ее границ — поверхности тела ОА ... Л7, изменения во времени параметров набегающего потока или параметров вдуваемого через поверхность тела на некотором участке АгА4 газа и т.д. Для тел сравнительно простой формы в возмущенной части течения вблизи затупления имеется единственная область дозвукового течения, за которой расположена область сверхзвукового течения. Если в некоторый момент времени t = 0 известно положение поверхности тела, ударной волны и распределение газодинамических функций, то течете в любой момент времени / > О полностью определяется заданием движения поверхности тела. Для определения этого течения необходимо найти решение уравнений (7), (9)-(11), удовлетворяющее начальным условиям при t = 0 и граничным условиям на поверхности тела и ударной волне, положение которой заранее неизвестно и подлежит определению вместе с течением.

Для более сложных течений, например, при сверхзвуковом обтекании затупленного тела, с части которого АъОА4 осуществляется интенсивный вдув газа (в общем случае другого состава, чем в набегающем потоке) в ударный слой, картина обтекания существенно усложняется. Между ударной волной 03Ci ...С5 и телом появляется поверхность контактного разрыва 02В - В7, отделяющая набегающий газ от вдуваемого и, так же как ударная волна, подлежащая определению. Основная область дозвукового течения А2В2В 3C2C3B4BSAS принимает сложную форму, за участком А3А4 вблизи тела могут возникнуть области возвратных течений, при очень сильном вдуве возможно появление сложной системы внутренних скачков уплотнения, поверхность контактного разрыва может замыкаться на поверхности тела и т.д.

Чтобы охватить по возможности все основные типы граничных условий, возникающих в задачах сверхзвукового обтекания тел, будем предполагать, что область возмущенного течения (ударный слой) состоит из области I, ограниченной контактным разрывом и ударной волной, области II между поверхностью тела и контактным разрывом, в обеих этих областях существует участок чисто сверхзвукового течения и что интенсивность других возможных разрывов существенно меньше интенсивности ударной волны и контактного разрыва.

Для определенности в дальнейшем также принято, что связанная с телом криволинейная ортогональная система координат (с началом координат — точкой 0 внутри тела) выбрана таким образом, чтобы в области интегрирования координата х2 возрастала в направлении потока в ударном слое, координата х3 соответствовала различным меридиональным плоскостям, а координата х2 имела направление, в целом соответствующее внешней нормали к поверхности тела (рис. 5.1).

Наиболее часто будет использоваться сферическая система координат х, =s, х2 =r, х3 = <р, связанная с декартовой х, у, z соотношениями

Рис. 5.1

В этом случае Яi - г, Н2 = 1, Н3 = г sin s.

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

  • 147
  • 173

В расчетах обтекания неосесимметричных тел при наличии угла скольжения р (рис. 5.1) координата 0<д:3<х^ = 2яи используется периодичность решения по этой координате

Удобно и в случае (13) использовать тот факт, что решение может быть известным образом продолжено за плоскости х3 = 0 и х3 = я, а именно

  • б) На ударной волне 03СХ ... С5 с некоторым уравнением х2 х3)
  • (которое, как отмечалось, подлежит определению) требуется выполнение законов сохранения массы, импульса и полной энергии, которые запишем в виде

— скорость набегающего потока газа относительно поверхности ударной волны в направлении ее нормали;

— наклоны ударной волны к осям х%9 х3 соответственно; орты выбранной системы координат.

Недостающие соотношения для определения входящих в (15) неизвестных р, Л, р, v2, в, qlc, q3c, Rc можно получить, привлекая уравнение состояния (7), три очевидных геометрических соотношения [44] (см. также разд. 2 гл. И)

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

явно выражается через в.

в) На участке поверхности тела А ... А3Л4 ... А7 с известным уравнением *2 = RwV, xi> *з) в отсутствие вдува граничным условием является условие непротекания

  • — единичная нормаль к поверхности тела.
  • г) На участке А2а поверхности тела со вдувом, в зависимости от физической природы появления вдува (истечение газов из сопла в ударный слой, разрушение теплозащитного покрытия под действием больших тепловых потоков и т.д.) возможны различные постановки граничных условий. При этом следует учитывать ограничения на допустимый вид граничных условий, следующие из разд. 5,гл. I. В частности, если механизм вдува обусловлен достаточно большой величиной теплового потока на единицу поверхности тела q(t, Х, дг3), вдув осуществляется по нормали nw к поверхности тела, т.е. V = Vnw = Vit^ + иэе^, и

является дозвуковым (V/а < 1, где а(р, h) - скорость звука), все термодинамические процессы равновесны и известно соответствующее уравнение Клайперона-Клаузиуса для зависимости давления насыщающих паров от температуры или энтальпии

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

(р, — плотность материала тела в твердой фазе), то такая простейшая модель для определения искомых параметров р, р, h, vlt v2, v3 (в случае, если D < V и Rw(t, Х, х3) — известная функция) включает уравнения (7), (22) и следующие соотношения:

Последнее выражение выписывается на основе анализа баланса энергии на поверхности тела (волне разрушения) и именно его различные формы отличают те или иные модели разрушения. В более сложных моделях часто учитывается также неравновесный характер процесса разрушения (в тонком пристеночном Кнудсеновском слое) и более детально - процессы в твердой фазе. Для замыкания системы (7), (22), (24), (25) следует привлечь определенную линейную комбинацию уравнений газовой динамики (9). В тех ситуациях, когда нарушается условие дозвукового характера вдува, в качестве замыкающего условия следут положить

так как ’’естественный” вдув не может быть сверхзвуковым.

Для определения формоизменения поверхности тела, величины уноса массы и т.д. необходимо привлекать закон сохранения массы (23) для Nw = D:

из которого можно определять функцию Rw(t, xit х3) и ее пространственные производные.

Вместо (22), (25) могут быть заданы величины

  • (в дозвуковом случае ’’принудительного” вдува) или даже полностью вектор и в случае сверхзвукового вдува
  • д) На контактном разрыве В если искомое уравнение его

поверхности есть х2 = Rk(t,xi,xз), нормаль nk = (-qikexl + ех2

- <7з*ех,)//1 + Як + Язк> наклоны к осям х х3, соответственно Чк = 2/Нi)(dRk(t, xif х3)/дх1), q3k = (Я23)(ЭДк/Эх3), газодинамические параметры по обе стороны от этого разрыва Uj и иц, то граничными условиями будут условие непротекания

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

е) При численном решении задачи об обтекании передней части затупленного тела область интегрирования ограничивают обычно некоторой замыкающей поверхностью {В iCiCsC^B 6А6А2). При определенном ее выборе (пространственный тип этой поверхности [17,18]) на ней не требуется постановка граничных условий. Условием пространственности этой поверхности является, например, выполнение в каждой ее точке в каждый момент времени неравенства

или, в координатной записи,

и стационарных параметрах вдуваемого газа, также можно использовать описанную выше постановку задачи. В этом случае в качестве начальных данных выбирают некоторые довольно произвольные функции R?t

Rc у u°, вообще говоря не удовлетворяющие (9), и отыскивают решение (9), (7) при стационарных условиях в набегающем потоке газа (37) и стационарных условиях вдува. Тогда в пределе при / -*00 возможен выход на искомое стационарное решение, т.е. отыскивается решение вида

Этот прием решения некоторых стационарных задач называют методом установления и широко используют в вычислительной практике.

3. Выбор системы координат в сравнительно простых случаях не имеет принципиального значения и определяется в основном соображениями удобства. Однако в ряде случаев область течения может иметь сложную геометрию со многими поверхностями разрыва, зонами больших градиентов и т.п. В этом случае либо используют сложные координатные системы, адаптирующиеся к картине течения и связанные с поверхностями разрыва- (например, [116, 117] и др.) или с искомым решением во всей области (например, [118, 119] и др.), либо строят однородные консервативные схемы [68], позволяющие с необходимой точностью сквозным образом рассчитывать течения при наличии особенностей внутри области интегрирования (см., например, [19]). По-видимому, наиболее эффективный путь в подобных ситуациях — сочетание обоих этих подходов с явным выделением основных особенностей течения (головной ударной волны и т.п.) и сквозным расчетом разрывов сравнительно небольшой интенсивности (внутренние скачки уплотнения, волны разрежения и т.д.). Не останавливаясь на детальном обсуждении этих вопросов, представляющих самостоятельный интерес, воспользуемся здесь простейшей заменой независимых переменных типа (1.2.10), переводящей криволинейную область интегрирования в более удобный для численных расчетов прямоугольный параллелепипед с фиксированными границами (при выбоце в качестве замыкающей повеохности ('ЗЗ'П:

Окончательно сформулируем рассматриваемую задачу. В области

найти решение (40), (7) и функции Лс(/,Я* (7 >?>?)> удовлетворяющие граничным условиям (13) или (14), (15)-(20), (21), (29) — (30), а также (22) —(25) или (24), (28) и некоторым начальным условиям (35),(36).

4. Вводя в области (42) равномерную по координатам ?,*?,? из (39) разностную сетку (рис. 5.2)

и обозначая значения искомых функций в ее узлах через ъпт1к> R”mk> Rпктк»••• Для расчета внутренних узловых точек, можно воспользоваться любой из рассмотренных ранее явных разностных схем, в частности, недивергентной схемой I из гл. III, т.е. (4.4.4). Для этого достаточно найти соответствующие собственные числа и собственные векторы матриц А|—Аз.

В частности, для матриц А х - Аэ из (9), (41) имеем

Непосредственной проверкой можно убедиться, что А/ =il~' AjQ./(i = = 1, 2, 3), где Л/ = ;Х/( - диагональные матрицы из собственных значений матриц А/.

Для построения расчетных формул на поверхности ударной волны

из которых формируются замыкающие условия для расчета граничных точек. Верхний знак берется при X2 >0, нижний - при Х;2 < 0.

При г? = 1 (узловые точки т = 0, 1,... ,М, / = L, к = 0,..., К), учитывая, что скорость распространения ударной волны по ’’жидким’* частицам в направлении своей нормали больше местной скорости звука за ударной волной, т.е. что -а < 0 <0, собственные значения матрицы Л2 имеют следующие знаки:

и в качестве недостающего соотношения к разностной аппроксимации граничных условий (15)-(19) используем (46) при 1 = 1. Разностная аппроксимация (19) служит для определения геометрических параметров поверхности ударной волны Уэстк* а подстановка (15) в

(7), (46) дает в общем случае нелинейную относительно систему из двух уравнений в каждом сеточном узле на поверхности ударной волны

которая решается обычным образом, т.е. методом итераций типа метода Ньютона.

На участке поверхности тела без вдува (rj = — 1, / = 0) с учетом условия непротекания (21), которое в новых переменных принимает вид и = 0, имеем

т.е. для расчета таких точек необходимо воспользоваться соотношениями (46) при I = 2, 3, 4, 5 и разностной аппроксимацией условия непротека- ния (21). В результате получим в каждом сеточном узле, принадлежащем этому участку поверхности тела, линейную относительно искомого вектора umo* систему из пяти уравнений, из которой несложно получить явные расчетные формулы для компонент этого вектора.

На участке поверхности тела со вдувом v = —H2R^t + V/l + q w + qw > > 0, следовательно, при дозвуковом вдуве

т.е. для расчета этих точек необходимо воспользоваться одним из соотношений (46) при / = 5 и разностной аппроксимацией граничных условий (22)-(25) или (24), (28), а при звуковом или сверхзвуковом вдуве - полностью задать вектор

В каждой точке контактного разрыва (т? = 0, / = L/2, т = 0, 1, . . . ,М, к = 0, 1, . . . , К) необходимо определить новое значение R^mk и Ю компонент векторов (и^д,/2,*)ь (um+,I/2,fc)n* ® каждой такой точке имеем v =vu = 0, т.е.

и для их расчета необходимо использовать 8 соотношений (46) при / = = 2, 3, 4, 5 в области I и при i = 1, 2, 3, 4 в области II, а также разностную аппроксимацию граничных условий (30), (31). Разностная аппроксимация (31) служит также для определения . Последующим численным дифференцированием R^k (или> аналогично, (19)) можно определить пространственные Производные Я*ктк>Чъктк-

Точки, принадлежащие замыкающей поверхности ЛХВХ СС^С^ВьАьАп

  • (т = М, I =1,2.....L - 1, к = 0, 1, . . . , К), при выполнении условия
  • (34) рассчитываются, как и внутренние узловые точки, поскольку X! = = (uj + а)/Н1 > 0, А/ = vx/Hx > 0 (i = 2, 3, 4), XJ = (у, - а)/Нх >0, и в соотношения (4.4.4) войдут значения и на слое t = tn = пт только в узловых точках 1 М - 1, /, к I, МУ I ± 1, к , i Л/, /, к ± 1 J, i М, /, А:) , т.е. в граничных и внутренних точках области интегрирования.

Расчет узловых точек к = 0 или к = К при m = 1,2,..., М - 1, / = 1,... . . . , L — 1, принадлежащих поверхностям х3 = 0 и х3 =xj, также производится по формулам для внутренних узловых точек с учетом граничных условий (13) или (14). При этом вводятся два дополнительных слоя узловых точек к = — 1 н к - К + 9 параметры в которых определяются по данным в симметричных точках к = ик=К-ъ соответствии с граничными условиями. В частности, для условия (13) имеем

В узловых точках, принадлежащих последней грани параллелепипеда (42) т = 0, / = 1, . . . , L — 1, к = 0, 1, . . . , К, при выборе системы координат (12) имеется особенность, связанная с этой системой координат. Для ее устранения при расчете этих точек делается формальный переход к сферической системе координат, повернутой относительно основной на я/2 вокруг оси у (рис. 5.1), в которой производится расчет узловых точек m = 0, / = 1, . . . , Z, — 1, А: = 0. Параметры в точках т = 0, / = 1,...,/. — 1, к = 1, . . . , К определяются с учетом того факта, что, по существу, они принадлежат одному и тому же физическому лучу 03 (jcj = 0).

Следует обратить внимание, что в нескольких ближайших по направлению Xi к этому лучу узловых точках m = 1,2,... 3)т1 к = х2/ sinxlm ^ **mhiX2i и для обеспечения аппроксимации исходных уравнений разностными выражениями (4.4.4) необходимо выполнение условия И3 h, что не требуется по реальному поведению решения в зависимости от координаты х3. Для устранения этого ограничения по направлению х3 в этих точках использовалась аппроксимация со вторым порядком точности, что обеспечивало в целом первый порядок точности, включая ближайшие к лучу *i = 0,лг3 =0 узловые точки.

 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ     След >