Лекция №6. Метод Распада Произвольного Разрыва (Линеаризованный)

Наиболее универсальными с точки зрения организации расчета нестационарного течения являются конечно-разностные методы, использующие фиксированную сетку в системе координат расстояние – время. К ним относится, в частности, метод распада произвольного разрыва (МРПР), основные положения которого разработаны С.К.Годуновым.

МРПР, как и другие численные методы газовой динамики, опирается на систему базовых интегральных уравнений, выражающих основные законы сохранения – массы, импульса (количества движения) и энергии:

; (215)

; (216)

. (217)

Здесь: х, t - система координат расстояние-время, р, ρ, Т, v, R, c_p, c_v, е = u+v²/2, λ _т, α _w, T_с, D - соответственно давление, плотность, температура, скорость потока, газовая постоянная, изобарная и изохорная теплоемкости, удельная энергия (сумма отнесенных к единице массы внутренней u = c_vТ, и кинетической v²/2 энергий), коэффициенты трения и теплоотдачи в стенку канала, температура стенки канала и гидравлический диаметр канала (в случае круглого трубопровода равен его внутреннему диаметру).

В данной постановке принимается физическая модель термодинамически идеального и совершенного газа (p = rRT, c_p = соnst, с_v = соnst, R = c_p - c_v = c_v(к-1), для воздуха c_p = 1005, 55 Дж/(кг·К), с_v = 718, 25 Дж/(кг·К), R = 287, 3 Дж/(кг·К), к = 1, 4) без воздействия внешних сил и источников тепла. Для учета диссипативных явлений используются представления Дар­си-Вейсбаха (вязкое трение) – правая часть уравнения (216) и Ньютона (теплоотдача в стенки канала) – правая часть уравнения (217).

При моделировании течения холодного воздуха (впускной тракт) теплообменом со стенками можно пренебречь, и правая часть уравнения (217) принимается равной нулю. В первом приближении расчетов можно также отказаться и от учета трения – принять равной нулю и правую часть уравнения (216).

Для организации расчета трубу длиной L разбивают на m отрезков. Тогда шаг расчета по координате х, т.е. линейный размер Dx расчетной ячейки составит

. (218)

Необходимо также выбрать шаг расчета по времени Dt. Расчет будет устойчивым, т.е. выполняться без сбоев, если при выборе Dt будет выполнено условие для так называемого сеточного числа Куранта Со:

, (219)
где а₀ – начальная скорость звука. Для нормальных условий с холодным воздухом:

= 347, 4 м/с.

Как правило, для обеспечения устойчивости расчета описываемой модификацией метода выбирают Со = 0, 5-0, 6. Тогда для шага по времени при расчете течения воздуха будем назначать:

. (220)

Теперь можно выстроить ортогональную расчетную сетку из m ячеек по горизонтали (по оси х) и N ячеек по вертикали (по оси y) (рис.41). Ячейки в горизонтальных рядах нумеруются: х₁, х₂,..., х_i,..., х_m. Границы ячеек нумеруются дробными индексами, например, для правой границы i – ячейки имеем х_i+1/2. Число ячеек по вертикали получится автоматически, если назначить окончание времени расчета t_N:

. (221)

Расчетные параметры снабжаются также и верхними индексами, соответствующими временным слоям сетки. Эти индексы меняются в интервале 0,..., n,..., N, тогда, например, если имеем , то это давление в i – ячейке в момент времени tⁿ.

Для организации расчета должны быть заданы начальные (НУ) и граничные условия (ГУ).

В качестве НУ необходимо назначить распределение всех расчетных параметров по всем ячейкам m ячейкам первого временного ряда в момент времени t⁰ = 0.

Например: v⁰₁=... v⁰_i...= v⁰_m = 0 м/с, р⁰₁=... р⁰_i...= р⁰_m = 100000 Па, Т⁰₁=... Т⁰_i...= Т⁰_m = 300 К, r = p/RT = ρ ⁰₁=... ρ ⁰_i...= ρ ⁰_m = 1, 160 кг/м³.

В качестве ГУ слева, т.е. в начале координат, будем рассматривать открытый конец со скруглением, которое обеспечивает втекание в трубу без потерь. Тогда, считая процесс адиабатным и не прибегая к ГДФ, можно воспользоваться уравнением энергии. И если в объеме слева имеем давление р_а = р_а* и температуру Т_а = Т_а*, для связи со скоростью и давлением втекания будет справедливо:

. (222)

Соответственно для температуры и плотности:

; (223)

. (224)

В качестве ГУ справа, т.е. в конце трубы, будем рассматривать либо открытое проходное сечение, когда давление на границе равно давлению в полости истечения, в простейшем случае – атмосферному:

, (225)
либо закрытый конец, когда скорость на границе полностью затормозится, т.е. станет равной нулю:

. (226)

Расчет выполняется последовательно от 1 до m ячейки. Сначала для всех ячеек определяются объемные параметры в исходный (начальный) момент времени t⁰ = 0: полная (т.е. не удельная) масса

, (227)

полное количество движения

, (228)

и полная энергия

. (229)

Если принять единичное значение площади F проходного сечения трубы, то в качестве объема ячейки можно принять ее линейный размер, т.е. шаг по длине

. (230)

В общем случае параметры v⁰_i, р⁰_i, Т⁰_i, ρ ⁰_i и, соответственно, объемные параметры М⁰_i, К⁰_i, Е⁰_i могут быть во всех ячейках ряда различными, т.е. будем иметь дискретное распределение параметров по длине трубы. Тогда труба может быть представлена в виде отдельных объемов Δ х, заполненных газом с разными параметрами, и которые по границам разделены перегородками (мембранами).

Таким образом, имеющееся на практике гладкое распределение параметров по длине заменяется разрывным, на каждой из границ имеется произвольный разрыв параметров.

В каждый из расчетных моментов времени tⁿ, в том числе и в начальный момент, когда tⁿ = t⁰ = 0, происходит как бы мгновенное разрушение всех мембран. В результате на границах между ячейками (дискретными объемами) происходят распады произвольных разрывов (РПР). На каждой из границ образуются две элементарных волны, фронты которых распространяются один —вдоль С⁺ ‑ характеристики, идущей от точки разрыва (границы) вправо вперед, другой —вдоль С^- ‑ характеристики, идущей от точки разрыва влево назад (рис.2).

Если, например, давление в i - ячейке в n - момент времени было больше давления в i+1 - ячейке, а скорости были равны нулю, то в этот момент времени от границы i+1/2 влево, внутрь ячейки i, пойдет волна разрежения, а вправо, внутрь ячейки i+1, - волна сжатия.

Рис.42. Образование элементарных волн при РПР

В соответствии с теорией нестационарного течения, при переходах через фронты элементарных волн выполняются законы сохранения инвариантов Римана. Так при переходе через фронт, идущий слева направо, т.е. вдоль С⁺ характеристики, будет оставаться постоянным инвариант Римана r^- с выполнением соотношения

. (231)

При переходе же через фронт, идущий справа налево, т.е. вдоль С^- —характеристики, будет оставаться постоянным инвариант Римана r⁺ с выполнением соотношения

. (232)

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

; (233)

, (234)

где

- (235)

осредненное значение произведения плотности и скорости звука для решения задачи РПР на i+1/2 границе.

Из системы уравнений (233)-(234) легко выразить «потоковые» давление и скорость через границу i+1/2 при РПР:

; (236)

. (237)

Для определения потоковой плотности используются зависимости

; (238)

, (239)

которые представляют собой линеаризацию адиабатных соотношений.

Для расчета на левой границе из (233) будем иметь:

. (240)

Поскольку с другой стороны имеется формула (222) для расчета скорости через границу, приравнивая (222) и (240), получим расчетное уравнение

. (241)

Решая его одним из методов последовательных приближений (практика показала, что надежнее всего использовать метод деления отрезка пополам), находим давление p_1/2 втекающего через границу газа. Затем по формуле (240) найдем соответствующую скорость, по формулам (223) и (224) – температуру и плотность.

Рассмотрим правую границу, где последняя ячейка расчетной области будет иметь индекс m, а узел непосредственно на границе - m+1/2. При РПР на этой границе налево вдоль С^- - характеристики пойдет элементарная простая волна, при переходе через фронт которой инвариант r⁺ остается постоянным, и из уравнения (234) можно получить в конечных разностях:

. (242)

Если конец трубы закрыт, то, очевидно, v_m+1/2 = 0. Это положение позволяет определить давление на границе:

, (235)

Затем:

. (236)

Здесь также можно отметить, что формула (235) отражает тот факт, что при подходе к закрытому концу волны сжатия (т.е. когда в данном случае v_m > 0) давление на границе в результате РПР возрастает, а при подходе волны разрежения (v_m < 0) упадет ниже первоначального р_mⁿ.

Пусть труба справа открыта, и происходит истечение. Тогда (225), и при подстановке в (234) получим скорость истечения:

, (237)

откуда видно, что при подходе волны сжатия (p_m > p_а) скорость истечения v_m+1/2, направленная по оси координат х, превзойдет скорость v_m.

Для плотности будем иметь:

. (238)

Обратимся теперь к исходным интегральным уравнениям (215)-(217). Решая задачу без диссипативных членов, т.е. без правых частей, отметим, что левые части представляют собой интегралы по контуру. И в конечных разностях для каждой из ячеек представляют собой сумму объемных параметров, полученную при обходе ячейки по ее прямоугольному контуру против часовой стрелки. Тогда можно записать:

; (239)

; (240)

. (241)

Эти уравнения показывают, что в трубе с единичной площадью проходного сечения объемные параметры: полная масса газа, полное количество движения и полная энергия объеме i ‑ ячейки за расчетный шаг по времени Dt, т.е. между tⁿ и tⁿ⁺¹ временным слоями численного расчета изменяются за счет потоковых параметров на левой (i‑ 1/2) и правой (i+1/2) границах ячейки. Т.е. при положительном значении скоростей через левую границу прибудет масса М_i-1/2 = (ρ v)_i-1/2Δ t, количество движения К_i-1/2 = (Мv)_i-1/2 и энергия Е_i-1/2 = (Ме)_i-1/2, через правую границу произойдут соответствующие вычитания М_i+1/2 = (ρ v)_i+1/2Δ t, К_i+1/2 = (Мv)_i+1/2 и Е_i+1/2 = (Ме)_i+1/2.

Кроме этих составляющих на изменение количества движения по сравнению с исходным значением окажет влияние положительное воздействие импульса силы давления на левой границе I_i-1/2 = p_i-1/2Δ t и отрицательное - воздействие импульса силы давления на правой границе I_i+1/2 = p_i+1/2Δ t. На изменение энергии по сравнению с исходным значением окажет влияние положительная работа, совершенная силой давления на левой границе L_i-1/2 = (pv)_i-1/2Δ t, и отрицательная работа, совершенная силой давления на правой границе L_i+1/2 = (pv)_i+1/2Δ t.

После определения объемных параметров М, К и Е в n+1 момент можно определить и основные газодинамические параметры в этот момент:

; (242)

; (243)

; (244)

. (245)

Очевидно, из последних двух формул можно выразить непосредственно:

. (246)

После завершения расчета по всем i – ячейкам ряда, необходимо перейти к следующему временному слою, заменив верхние индексы n+1 всех параметров на n, и так далее, переходя к последующим рядам вплоть до окончания назначенного расчетного времени t^N.

Расчет временного ряда начинается с определения потоковых параметров при входе через левую границу (241), (240). Далее следует расчет внутри расчетной области по формулам (236), (237), (238) или (239) по всем i, наконец, правая граница (235), (236), если конец закрыт, и (237), (238), если конец трубы открыт. Затем определение объемных параметров:

; (247)

; (248)

, (249)

где . (250)

Расчет может быть организован вообще без использования температуры. Поскольку ,

. (251)

Затем следует повторение балансовых расчетов (239)-(241).