К выбору метода численного интегрирования уравнений одномерного движения газов в каналах ДВС

№6

2007

621.436

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

В КАНАЛАХ ДВС

Канд. техн. неук, доц. А.А. ЧЕРНОУСОЕ

Три разностные схемы для решения уравнений одномерной нестационарной газовой динамики сравниваются по виду численных решений тестовых задач о распаде разрыва и по затратам машинного времени. Первая — классическая схема С. К. Годунова, две другие — схемы того же типа повышенного порядка аппроксимации. По вычислительной эффективности выделяется одношаговая схема повышенного порядка аппроксимации по пространству и второго порядка по бремени. Схема рекомендована к применению для расчетов течений в каналах газовоздушных трактов ДВС.

Three finite-difference schemes for solving one-dimensional conservation laws of gas dynamics are compared by numerical solutions obtained and by computational effort required. The first one is the classical Godunov scheme, other two are higher-order Godunov-type schemes. One-step scheme of second order of accuracy in time is found to be the best in terms of numerical efficiency. This scheme is recommended for applications such as numerical modeling of unsteady gas motion in channels of internal combustion engines.

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

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

Имеется много методов численного интегрирования систем квазилинейных гиперболических уравнений. Среди них выделяются явные консервативные монотонные методы сквозного счета повышенной точности, позволяющие получать численные решения с высоким разрешением и без нефизичных осцилляций. Задача выбора наилучшего метода данного класса возникла в связи с модернизацией расчетного модуля ТРУБКА системы имитационного моделирования «Альбея» [1]. Необходимо было выбрать метод интегрирования уравнений, обладающий, наряду с консервативностью и монотонностью, также высокими порядком аппроксимации и вычислительной э ффективн о с т ы о.

Такой метод отыскивался в подклассе методов типа С.К. Годунова, к которому условно относят разновидности метода конечных объемов, где газодинамические потоки на границах ячеек определяются из решения локальной задачи о распаде разрыва точного или линеаризованного вида или из локальной характеристической задачи.

№ б

2007

1. Законы сохранения и разностные схемы 1.1. Законы сохранения одномерной газовой динамики

В основе описания неустановившегося движения газовой смеси в канале как квазиодномерного лежит система интегральных законов сохранения массы каждого из К компонентов, количества движения и полной внутренней энергии смеси, записанных для конечного участка канала Ах

-1 V.

А* Г

Э/ J

рkAdx = -(pkuA) k - U.... К

Э/

J р и A dx

р if + р

и

Л-., -1 Дг

A'n -i Vv

+

ГТ dA

Хп.\\ +[j — dx

dx.

0_ Э/

, ! Av

л,, I .U

J рЕА dx - - [(pw£ + ри) А ] | " + J qtvY\ dx

где искомые функции— рДх, t), к = \,..., К, и(х, /) и Е(х, /), А-А(х) — плошадь сечения канала, П = П(д-) — его периметр, т„, и qlt, — средние но периметру величины касательных напряжений и тепловых потоков на стенке.

Температура и давление связываются с искомыми функциями уравнений (I) — (3) через уравнения состояния (УС). Для смеси идеальных газов с помощью калорического УС температура Т(х, 1) находится по массовому составу Yk = рк /р, к = К и удельной внутренней энергии е(х, /) смеси из равенства Е = ^Г , ekQ Wk +0,5if , а давление

р(х, г) —из термического УС р-рЯТ, где Л" =

В подобластях гладкости (дифференцируемость по независимым переменным х и /) искомых функций р4. и, Т, р системе (1) — (3) эквивалентны уравнения в дифференциальной форме, получаемые из них в пределе Ах —» 0 :

+ = Лг = 1.....if,

at ох

д 1

д (риА) + —- Г(рк2 + р)А1 = хиЛ + р

ах

~ (р ЕА) + — [fp иЕ + ри)А~] = qw П, at ахL

dA dx

или, в оолее компактной записи, с применением символических «векторных» обозначений

где U = [р,,.... рл,, ри, р£']Г и F = [p,w,..., рКи, ри2 + р, риЕ + fti<~j — соответс-

твенно «векторы» неизвестных и плотностей газодинамических потоков, а

№ б 2007

- г <*А г/

5 = [0,О, хуу1\ +р — ,днЛ

— «вектор» объемных плотностей источников масс ком-

понентов, количества движения и энергии.

Форма записи уравнений (4) — (6) именуется дивергентной (или консервативной), чем подчеркивается, что в ней сохранено унаследованное от уравнений в интегральной форме представление о сохранении массы, импульса и энергии при изменениях их количеств в контрольных объемах под действием потоков и источников. Векторные обозначения вида (7) удобны и применяются ниже для записи соотношений численных методов.

При применении уравнений (4) — (6) или (7) для вывода методов численного интегрирования получаются консервативные разностные методы сквозного счста так же, как и при использовании для этого исходных уравнений (1) — (3).

1.2. Метод «распада разрыва» С.К. Годунова

За основу при сравнении погрешностей численных решений и затрат времени на счет взят классический метод «распада разрыва», или конечно-разностная схема С.К. Годунова [2, 3].

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

ин+х =£/."+- А/

АхА.

+ д. '

где - и Р.'^ = Р {у", ) — «векторы» газодинамических потоков из

решения задачи о РПР на границах /-ой ячейки (с индексами и /' + {) на «старом» временном слое. Данная схема (упоминаемая ниже как С1) была реализована в программе для ЭВМ с использованием процедуры решения задачи о РПР итерационным методом Ньютона из [3].

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

1.3. Схема типа Годунова повышенной точности

Недостатками классической схемы Годунова являются ее низкий порядок аппроксимации (и следствие — сильная схемная диссипация) и относи тельно булыние вычислительные затраты, а достоинством — монотонность решений. При простом расширении сеточного шаблона в целях повышения точности наблюдаются нефизичные дефекты решений, главным образом, осцилляции сеточных функций {немонотонность).

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

№ б

2007

-<:>........

Рис. 1. Схема обновления решения в ячейках по методу «распада разрыва» С.К. Годунова (С1); профили

решения типичны для плотности р

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

Схема повышенного порядка аппроксимации С2 строится на квадратичном представлении решения по ячейке (рис. 2). Реконструкция значений параметров на границе для решения локальной задачи о PHP использует процедуру диагонализации матриц системы гиперболических уравнений. Процедура расчета U по обе стороны, например, от / + \ -ой границы, имеет следующий вид:

= и;+[s-1 ]'; a w* + ^ V Щ" |, = й°м-[5-'];;, j^v^i.+^Äcj,

Щ" = min mod(AW", bVW"), VIV," = min mod(VIV,", bAlV/'),

V W" = [S]" V Ö" - [>S]" (Ü¡ - ¿7,",),

где 1 < b < ömax, b = -——, ф < 1, а ограничительная функция min mod(-, •) определяется как * mJX 1-ф

№6

2007

mm

in mod

[ 0, ху < О

[sign(x) min(| х\,\у\), ху> О

Значения параметров процедуры, обеспечивающие третий порядок пространственной аппроксимации в аналогичной схеме для модельного линейного гиперболического уравнения: (р = 1/3 , Ъ = Ьтт = 4 ([6]).

MiiiStS

iiliiii^eire

I г

Рис. 2. Обновление решения ячейке по схеме повышенной точности С2 (двухшаговый вариант «предиктор-

корректор»)

Матрица [5] и обратная ейматрица [5-1] взяты из преобразования приращений вектора консервативных переменных и в приращения вектора плотностей потоков ^ , справедливого для системы уравнений одномерного плоского движения однородного совершенного газа или смеси совершенных газов постоянного состава (у = с /су):

= М] = [5-,][Л][5], [Л] = с11а%{и,и + с,и-с},

№ 6

2007

[5]:

> у-1 о

- с" +-гг (1 - у)и у -1

У-1 о

- си н--ы с + (1 - у )и у -1

у -1 ,

С14 +-и" -С + (\-у)и у-1

[<Г]] =

с и

2с

и 1

—т +— 2с2 2с

2с и 1

и" и - +— +

2с2 2с и

--+ -

2с2 4с2 2с 2

Для расчета потоков массы компонентов смсси на границе необходимо интерполировать на нее набор парциальных плотностей рА.. Это легко сделать, дополнительно применив процедуру описанного выше вида для массовых долей Ук; для них, как характеристических переменных системы (4) —- (б), матричное преобразование в процедуре реконструкции не требуется.

Второй или третий порядок аппроксимации по времени в схеме С2 достигается соответственно в двух- и трехэтапной схеме обновления решения па шаге по /. Двухэтатшая схема «предиктор-корректор» записывается как

2 2 2

а трехэтапная схема обновления решения в ячейке (метод Рунгс—Кутта), рекомендуемая В [7],

й1]) = и;' + А Ш'Фч),

' 4 4 4

и"+] = - и" + - и\2) + - At.ll1 (Р\2)). 3 ' 3 3 '

В этих схемах 1!'(У") — разностный пространственный оператор; применительно к (7) он имеет вид

АхА;

В описываемой схеме для расчета потоков по данным, интерполированньш на границы, могут применяться как процедура точного расчета РПР, так и приближенные процедуры. Более экономичная характеристическая процедура может быть построена на системе скалярных уравнений (производные берутся вдоль характеристических кривых) = 0, ¿/°/() = 0, ¿Г/, = 0 для инвариантов 1_ = и + а_р9 /0 = р + а()р, /,=г/ + а.р

№ б 2007

линеаризованной системы уравнений плоской одномерной газовой динамики в характеристической форме

d~u _ d-p = о, d[)p - с V V = 0, w + — d + р = 0,

рс рс

в которой а_ = -1 /(рс), а0 = -с2 и а, = 1/(рс) — постоянные коэффициенты.

Поясним способ вычисления параметров в некоторой точке D (например, на границе ячеек) по значениям инвариантов на прибывающих в точку характеристиках системы (8). Для определенности возьмем случай (г/ ~ с) < 0, и > 0, показанный на рис. 3:

Ld = ив + а_рв, /()D = рс + а0рс, /+ 0 = + а_рА.

Простейший способ расчета коэффициентов в инвариантах — а_ = а_в, а() = аос, и ol = а+/1; его и предпочтительно применять для дозвуковых течений.

Наконец, по известным инвариантам и коэффициентам а_, а0 и а+ вычисляются параметры в точке D

/+ - /_ , h-P

Р = —-, и = 1+-а+р, р = --?

а+-а_ а0

после чего по УС определяется температура Т = p/(pR), где постоянная i? или массовые доли Yj,YK для ее расчета берутся для той смеси, которая, согласно знаку скорости потока, течет через границу

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

1.4. Экономичная одношаговая высокоточная схема

Инварианты линеаризованной системы уравнений одномерной газодинамики дают возможность построить очень экономичную монотонную схему для численного решения одномерных нестационарных задач. Идея и описание такого класса схем содержится в [8]. Опишем такую схему для уравнений (7), обозначив ее СЗ.

В СЗ применяется процедура квадратичной реконструкции решения в ячейке без матричного преобразования _ _ _ __

7«/ Ч 7« A//7 + V/.V ч А/."-V/"

//<(х) = //< + ' | (х_Х/) + ф I i.{х_х )-

2 Ах Ах'

AI/' = min mod(A7/\ bVT"), V//7 = min mod(V//, ЬЫ"\

где в число зависимых переменных/ = [/_,/0,/_ Д, YK^]r включены и массовые доли КЛ компонентов газовой смеси.

Параметры смеси и газодинамические потоки на границах (рис. 4) определяются в СЗ с помощью характеристической процедуры, подобной описанной в разделе 1.3 и в [8]. В данном же случае для расчета потоков на границах z-ой ячейки нужно определить параметры потока в узлах сетки (/ ±±9 п + как в точке D на рис. 3, т.е. по значениям инвариантов и коэффициентов а в точках Л, В и С, определяя эти величины по кусочно-параболическим распределениям /."(*) —на «старом» временнум слое (рис. 4).

№ 6 2007

О

Рис. 3. Шаблон линеаризованной локальной характеристической задачи при (и -с) < 0, и > О

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

и"+] = и" +

м

АхА,

А/

А. *'

в которой потоки «на полушаге» /^'¿-вычисляются описанным выше способом, а источ-никовый член может быть вычислен (для сохранения второго порядка по времени) как

"/А1?, л,'->,

Т"

-/"У—

+.....ф——........-О.............О I............... О-......О..........0~;.......О—

--О..............О.....»

Рис. 4. Обновление решения экономичной одношаговой схемой

№6

2007

2. Результаты тестовых расчетов

С применением всех трех схем решены две тестовые задачи о распаде разрыва; выбор задач обусловлен необходимостью проверки возможности работы данных монотонных схем как методов сквозного счета. В условиях обеих задач задавалась нулевая начальная скорость однородного газа с у = ср Icv = 1,4 = const, а отношение как плотностей, так и давлений по обе стороны от начального разрыва составляло 2 и 20 для первой и второй задач соответственно. В решении первой задачи весь газ в возмущенной при распаде начального разрыва зоне движется с дозвуковой скоростью, во второй задаче развивается сверхзвуковое течение с М ~ 1,36 в волне разрежения.

На рис. 5 показано решение первой тестовой задачи схемой Годунова С1, а на рис. 6 — решение второй задачи. Расчетная сетка содержала 100 ячеек, сделано 100 шагов во времени с числом Куранта, равным 0,5. Численные решения на рисунках наложены на точные решения тестовых задач.

^чГТ

1.8 и

1.4

1.2

V ; О *

П "Л

* п

%

; ui

Г: \

тлхтт

20

40 60

Ячейки

ВО

100

Рис. 5. «Дозвуковой» распад разрыва, рассчитанный схемой Годунова

В решениях по схеме С1 заметно сильное «размазывание» решения, хотя положение разрывов, как и ожидалось, воспроизводится правильно. В решении второй, «сверхзвуковой», задачи, виден дефект, присущий всем схемам типа Годунова первого порядка — не-физичный скачок параметров на звуковой линии.

Схема С2 была применена в «трехшаговом» варианте, с линеаризованной процедурой решения задачи о РПР на границах ячеек. Полученные по С2 решения задач показаны на рис. 7 и 8. Видно, что при отсутствии нефизичньтх осцилляций численное «размазывание» значительно меньше, чем в С1, а положение разрывов, как и в С1, воспроизводится правильно.

Решения тех же тестовых задач экономичной одношаговой схемой СЗ показаны на рис. 9 и 10. Данная схема показала наибольшую «разрешающую способность». В решении

№ б

2007

заметна некоторая «шероховатость», которая, очевидно, является следствием примененных в СЗ линеаризаций; впрочем, этот дефект существенно проявляется лишь в решении «сверхзвуковой» задачи (рис. 10).

20 ; 18 ¡-15 14

12 г

10

2 Ь

20

60

80

100

Рис. 6. «Сверхзвуковой» распад разрыва, рассчитанный схемой Годунова

. 5. П. 515.П. I.и .;...!.£]. и П и • ¿.11и.:.Ш ; ! Ш.: ! }••:

20 40 60 80 ¿00

Рис. 7. «Дозвуковой» распад разрыва; схема повышенной точности С2

№6

2007

20 :

18

14

12

ч

20

40 60

Ячейки

ВО

100

Рис. 8. «Сверхзвуковой» распад разрыва; схема повышенной точности С2

2

1.8 !

Кб Ь

1.4

1.2 г

t v

20

40

60

80

100

'Н.'ПК!:

Рис. 9. «Дозвуковой» распад разрыва; экономичная одношаговэя схема СЗ

Три рассмотренные разностные схемы требуют (по порядку величины) одного и того же времени счета на ЭВМ. Так, затраты времени ЦПУ в задаче о «дозвуковом» РПР на сетке из 1000 узлов по пространству и по времени на ПК с процессором Intel Pentium III

на тактовой частоте 700 МГц по этим схемам составили; для С1 — 3,58 с, для С2_8,33 с

и для экономичной одношаговой схемы СЗ — 2,23 с.

№6

2007

20 ;

18

16

14

12

И) :

6 Г

20

60

80

100

Рис. 10. «Сверхзвуковой» распад разрыва; экономичная одношаговая схема СЗ

Выводы

Из трех рассмотренных разностных схем для одномерных уравнений газовой динамики одношаговая экономичная схема СЗ обладает наименьшей схемной диссипацией, судя по виду численных решений в окрестностях разрывов. Она же оказалось самой экономичной по затратам машинного времени, потребовав 26% от времени счета по трехшаговой схеме повышенной точности С2 и 62% от времени счета по классической схеме Годунова С1. С учетом высокой «разрешающей способности» схемы СЗ по оценке, выполняемой решением задач о РПР на разных сетках, СЗ эффективнее схемы О в вычислительном плане примерно на 2 порядка.

Некоторая «шероховатость» численных решений тестовой задачи, полученных схемой СЗ, позволяет ее рекомендовать лишь в приложениях, в которых этот дефект не проявится заметно. Так, например, волновое движение в газовоздушных трактах ДВС как раз отличается умеренной интенсивностью волн, и для таких расчетов схема СЗ может быть рекомендована.

Практика эксплуатации модернизированного расчетного модуля ТРУБКА системы имитационного моделирования ДВС «Альбея» [ 1 ] показала, что результаты моделирования процессов газообмена ДВС с применением схем С2 и СЗ различаются незначительно.

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

1. Система имитационного моделирования «Альбея» (ядро). Руководство пользователя. Руководство программиста /

B. Г. Горбачев, С. А. Загайко, И. В. Рудая, и др. //Учебное пособие. — Уфа: УГАТУ, 1995. — 112 с.

2. Годунов С. К. Разностный метод численного расчета разрывных решений гидродинамики. / Матем. сб. — 1959. — Т. 47. — Вып. 3. — С. 276—306.

3. Численное решение многомерных задач газовой динамики / С. К. Годунов, А. В. Забродин и др.. — М.: Наука,

1976. —400 с.

4. К о л г а н В. П. Применение принципа минимальных значений производной к nocí роению конечно-разностных схем для расчета разрывных решений газовой динамики // Ученые записки ЦАГИ. — 1972. — Т. 3. —М> 6. —

C. 68—77.

№ 6 2007

5. Иванов М. Я., Ни гм ату л лин Р. 3. Неявная схема повышенной точности для численного интегрирования уравнений Эйлера. //Ж. вычисл. матем. и матем. физ. — 1987. — Т. 27. —№11. С. 1725 1735.

6. Ива н о в М. Я., К р у п а В. Г., Н и г м а т у л л и н Р. 3. Неявная схема С.К. Годунова повышенной точности для интегрирования уравнений Навье-Стокса. // Ж. вычисл. матем. и матем. физ. — 1989. — Т. 29. — № 6. — С. 888—901.

7. Gottlieb S., S h и С. - W. Total variation diminishing Runge-Kutta schemes. // NASA CR-201591 1CASE Report

№ 96-50, ICASE, 1996. — Pp. 20.

8. Копченов В.ИДрайко A.H. Монотонная разностная схема второго порядка для гиперболических систем с двумя независимыми переменными. // Ж. вычисл. матем. и матем. физ. — 1983. — Т. 23. — № 4. — С. 848—859.

629.424:621.436+004.94

ПРИМЕНЕНИЕ МЕТОДОВ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАБОЧИХ ПРОЦЕССОВ ДИЗЕЛЯ ПРИ ИНТЕРПРЕТАЦИИ РЕЗУЛЬТАТОВ ДИАГНОСТИЧЕСКОГО ЭКСПЕРИМЕНТА

Д-р техн. наук, проф. В. А. ЛАШКО, канд. техн. наук AM. КОНЬКОВ

Предложен метод идентификации неисправностей дизеля с использованием имитационного моделирования рабочего процесса дизеля. Приведены результаты экспериментальной проверки метода, полученные для двигателя 4VD 12,5/9 (44 9/12,5) дизель-генераторной установки рефрижераторного подвижного состава

Method of fault identification in a diesel engine using model that accurately simulates the diesel combustion process is offered. Results of experimental inspection of this method, obtained for the drive 4VD 12,5/9 (4ch 9/12,5) of a diesel plant at the refrigerator rolling stock, are obsewed.

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