CC BY
71
Уфа : УГАТУ, 2009 &.....& <£/ ~ ' Т. 12, № 1(30). С. 197-210
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
УДК 621.43:519.8
А. А. ЧЕРНОУСОЕ
О ДОСТОВЕРНОСТИ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ ДВИЖЕНИЯ ВОЛН КОНЕЧНОЙ АМПЛИТУДЫ В ДЛИННОМ НЕРАЗВЛЕТВЛЕННОМ ТРУБОПРОВОДЕ С МЕСТНЫМИ СОПРОТИВЛЕНИЯМИ В ОДНОМЕРНОМ ПРИБЛИЖЕНИИ
Исследовано прохождение волн конечной амплитуды (ВКА) по длинному трубопроводу с местными сопротивлениями (МС) в виде простых шайб. Расчеты в одномерном приближении выполнены с привлечением зависимостей для потерь полного давления на МС вида о(Мт), полученных продувкой МС стационарным потоком. Расчетные зависимости p(t) в сечениях трубопровода сопоставлены с экспериментальными. Показана высокая степень их соответствия, чем подтверждена достоверность одномерного приближения для моделирования движения длинных ВКА по неразветвленному трубопроводу c МС. Нестационарные течения ; модели местных сопротивлений ; численное моделирование
ВВЕДЕНИЕ
При моделировании нестационарных течений сжимаемой жидкости в трубопроводах в одномерном приближении численно решаются уравнения механики жидкости и газа, записанные для квазиодномерных (для гладких участков - каналов) и «нульмерных» (для емкостей) распределений зависимых переменных, с привлечением моделей местных сопротивлений (МС) в нестационарном потоке, задающих условия совместности в узлах стыков указанных элементов трубопровода.
Математическая модель течения через МС позволяет определить параметры потока в узловых сечениях на временном шаге такого расчета. Классические модели МС [1] объединяют в себе характерные для одномерной газовой динамики квазистационарные соотношения непосредственно на стыках элементов трубопровода и соотношения на элементарных волнах.
В [2, 3] показана возможность довольно точного расчета в одномерном приближении волнового процесса в протяженном трубопроводе как без существенных МС, так и с МС, статические характеристики которого в потоке сжимаемого газа определяются по приближенной расчетной методике. Для получения стабильно достоверных результатов моделирования в этом приближении характеристики МС должны формироваться (в ответственных случаях) по результатам рационально организован-
ных продувок МС.
Для рассматриваемых здесь относительно простых («однопоточных») МС, располагаемых на стыке двух каналов или на стыке канала и емкости, модель течения непосредственно на МС обычно строится на допущениях о квазистационарности (О = со^^) и энергетической изолированности (ОН = соп81;2) процесса перетекания, потери полного давления между входным и выходным сечениями МС учитываются коэффициентами £ или о. Критериальное уравнение, выражающее потери полного давления для течения определенной сжимаемой среды через заданное МС в указанной постановке должно иметь, например, вид
а = а(Мт, ЯеТ), (1)
где а = р*ых /р*х £ 1 - коэффициент восстановления полного давления, Мт и Яет - числа Маха и Рейнольдса в контрольном сечении, относящемся к сечению примыкающего канала («трубы»), по возможности - на входе (МТ = Мвх, Яет = Яевх) в МС. При течениях газов с большими скоростями (и числами Яет) зависимость (1) «вырождается» до
а = а(Мт) (2)
и важным становится адекватный учет эффектов сжимаемости и возможного «запирания» в критическом сечении на МС, что существенно, например, при моделировании течения через ор-
ганы газообмена ДВС.
Цель представленного исследования состояла в изучении степени достоверности результатов, получаемых численными расчетами нестационарных течений газа в протяженных трубопроводах с «однопоточными» МС по методике, включающей применение рациональной модели МС в нестационарном потоке и характеристик МС вида (2), полученных по результатам статических продувок.
1. УРАВНЕНИЯ КВАЗИОДНОМЕРНОГО ТЕЧЕНИЯ И ИХ СЛЕДСТВИЯ
Система уравнений квазиодномерного движения газовой смеси, лежащая в основе описания неустановившихся течений рабочего тела тепловых двигателей в одномерном приближении, имеет смысл системы законов сохранения массы каждого из К компонентов, количества движения и полной внутренней энергии смеси, записанных для участка канала Ах.
Для численного расчета описываемых этими уравнениями течений на участках каналов в модели элемента (МЭ) ТРУБКА системы имитационного моделирования «Альбея» [4] используются высокоточные явные методы (см., например, [5]). Такие методы аппроксимируют на конечно-объемной сетке уравнения в интегральной форме или же их следствия - уравнения в дифференциальной форме, получаемые в пределе Ах ^ 0:
ЭРк^ + Эр^
- + -
Эt Эх dpuF Э(ри2 + p)F
= 0, k = 1, ..., K,
(3)
+
Эt Эх
ЭpEF Э(ри£ + pu) F
= twn + pdF, (4) dx
Эt
+
Эх
= q^n.
(5)
Здесь рк(х, 0, к = 1, ..., К - парциальные плотности компонентов, и(х ,0 - скорость, Е(х, 0 - полная удельная энергия, р(х, ¿) - термодинамическое давление. Для замыкания системы уравнений (3) - (5) необходимо задать площадь ^ = ^(х) и П = П(х) периметр каждого сечения канала, термическое р(р, Т, У1, ...) и калорическое е(р, Т, Уъ ...) уравнения состояния многокомпонентной смеси (Ук = рк/р = рк/^рк), а также указать способ вычисления и -средних по периметру сечений локальных значений касательного напряжения и теплового потока на стенке.
Система уравнений (3)-(5) относится к классу систем квазилинейных гиперболических уравнений в частных производных [6, 7] и мо-
жет быть приведена к так называемой характеристической форме, полезной при создании как эффективных методов численного решения на участках каналов, так и моделей взаимодействия нестационарного потока в трубопроводе с местными сопротивлениями.
Для частного случая однородной (К = 1) среды с нормальной сжимаемостью характеристическая форма уравнений системы (3)-(5) принимает следующий вид [6]:
1 d+р + d0 p 2 d0 p
= b
- c
pc
где
1 (Эр
b* = +-b + b +—I ^ lb, b0 = | J b
p 1 2 pc У Э^1 3'
= b0, (6)
Эр
где 5 - удельная энтропия, с - скорость звука и, в свою очередь,
b =-PUdF, b2 = ЬП
F dx pF
1 dq b3 = —-■ TdT
(7)
где dq = dqвн+ dqвнеш - удельная теплота, подводимая к частице среды как от внутренних источников (трение о стенку, скачки уплотнения), так и от внешних (теплоотдача в стенку).
Производные от искомых параметров среды и, р, р в (6) берутся на плоскости (х^) в характеристических направлениях, определяемых как
d х
= и(х, t) + c(х, t),
= u(х, t).
Как известно [1, 6, 7], семейства интегральных кривых уравнений (7) - траектории индивидуальных частиц среды и элементарных упругих возмущений, распространяющихся со скоростью звука относительно частиц.
В теории важен частный случай одномерного нестационарного течения, когда сечение канала постоянно (F = const), напряжения трения Tw и тепловые потоки на стенке qw нулевые, в изучаемой области отсутствуют сильные разрывы искомых функций, а начальное распределение удельной энтропии ,t0) однородно. Тогда правые части уравнений (6) обращаются в нуль, и однородным будет распределение ,^(х, t > t0) энтропии в изучаемой области (изоэнтропное течение), при этом соотношения вдоль характеристических направлений примут особенно простой вид. В частном случае идеального совершенного газа (р = pRT, e = cvT, где cv = cp-R = const;, у = cp / cv = const2) эти соотношения будут выражать постоянство инвариантов Ри-мана I± = 2cl(y-l) ± и вдоль соответствующих
u
0
х
характеристических направлении.
Изменение одного из инвариантов Римана I± в течение указанного вида на постоянную величину связно с движением по газу простой изо-энтропноИ волны конечной амплитуды (ВКА), несущей возмущение I+ или I. Изменения параметров газа на фронтах таких волн выражаются газодинамическими функциями (ГДФ) [1]. Эти ГДФ удобны для записи моделей взаимодействия с МС нестационарного потока в трубопроводе.
Так, отношение статических параметров к параметрам газа, заторможенного во фронте волны, несущей возмущение инварианта Римана I_ и движущейся по невозмущенному газу во, определяется набором ГДФ от М = u/c и у:
1
-,=a'(Mg) = -c 1 -g1 M
T=t (M, g)=[a(M, g)] 2 = ( 1 )2 T (i -g1 m )
, (8)
2g lg-1
p=p(m ,g)=[a(M ,g)]g-1 =-p
i
1 -g1 M)
2g '
g-1
При торможении одномерного потока волной, несущей возмущение I- и движущейся по нему влево, используется другая группа ГДФ:
c=a(M g)=T7g^>
1 + g-i M
T = t(M,g) = [a(M,g)] 2 = --1-, (9)
T (l + ^ M j2
P=p(M ,g)=[a (m ,r)f =-
p (l+1-1 M,
2g ■ g-1
Благодаря обратимости процесса изоэн-тропного торможения, данные ГДФ удобны для выражения отношений одноименных параметров состояния газа по обе стороны от фронта простой изоэнтропной ВКА; например, на фронте такой волны, движущейся по газу влево, справедливы выражения
__(M 2,g)
p(Ml, g)' Ti t(Mi,g)
p = P(M2, g) T
Pi
и т. д.
В заключение приведем сводку соотношений, используемых в моделях МС для описании стационарного квазиодномерного течения. Так, исключая производные по ^ из (3) - (5) получим, что в стационарном потоке однородной (К = 1) сжимаемой жидкости параметры в сечениях 1 и 2 канала связаны соотношениями
(puF )2 = G2 = Gi = (puF)i'
[(pu 2 + p)F ] 2 = J2 =
=Ji+J
dF + p —
dx
dx,
[(puE + pu)F] 2 = (Gh*)2 =
= (Gh*)i + J XX2 qJPdx.
(i0)
(ii)
(i2)
При часто принимаемом допущении об энергетической изолированности (дм, = 0) потока на участке между сечениями 1 и 2 уравнения (10)-(12) эквивалентны соотношениям
G2 = Gi' J2 = Ji + J
dF + p —
dx
dx,
(i3)
h* = hi*.
Для идеального совершенного газа (h = cpT, cp = const) третье соотношение (13) дает Т2 = T*.
Среди соотношений (13), связывающих параметры потока в сечениях 1 и 2, второе нетривиально. Поэтому для расчета течения через МС по (13) требуется вычислить интеграл по контуру МС, что затруднительно. На практике для расчета потока через МС вместо уравнения сохранения количества движения используют связь полных давлений в сечениях, получаемую в общем случае продувкой МС и задаваемую критериальным уравнением, например, вида
p* = Si2(Mi, Rei, ...). pi
(i4)
Полное торможение стационарного потока совершенного газа теоретически может быть выполнено в обратимом изоэнтропном (йд=0, о=1) процессе. Отношения статических параметров в сечении к параметрам, получаемым в результате такого торможения, выражаются известными ГДФ [1, 9]
T i
T = t(M,g) = i+^M '
-cr = a(M, g) = [t(M, g)]2 = -
(i + ^ M2)
(i5)
4- = р(М, у) = [т(М, у)Г =-г.
Р (1 + ^ М2 )у-1
Данная группа ГДФ служит для выражения связи одноименных параметров газа в разных сечениях. Для энергетически изолированного стационарного потока совершенного газа
x
2
x
x
2
x
i
c
(Т2 = Т1) с учетом потерь полного давления (р2 = &12.Р 1) на участке 1-2:
Р2 Р\
= &
2, у) = Т(М2, У) я(Мх, у), Т т(Мх, у)
и т. п.
2. МОДЕЛИ МЕСТНЫХ СОПРОТИВЛЕНИЙ В НЕСТАЦИОНАРНОМ ПОТОКЕ
Для корректного расчета неустановившегося течения в трубопроводе с МС квазистационарных соотношений в модели МС недостаточно. Уравнения (3)-(5), используемые для описания течения на «гладких» участках каналов, образуют систему уравнений в частных производных гиперболического типа [6], что налагает на модель МС требование учета (хотя бы с линеаризацией) соотношений на прибывающих к сечению МС характеристических кривых в (х , ¿).
Этим требованиям отвечает, например, определение потоков через МС на расчетном шаге из решения обобщенной задачи о распаде произвольного разрыва (РПР) параметров в исходных данных из примыкающих к МС каналов и емкостей. Данный подход впервые применен, по-видимому, в [8], где задача о РПР на скачке сечения канала «замкнута» с привлечением гипотезы о среднем давлении на стенке в уравнении сохранения количества движения, записанном для контура, охватывающего скачок сечения трубопровода, а для волн сжатия, возникающих при РПР, применены соотношения на скачке уплотнения.
Изложенные далее модели МС основаны на близком подходе, развитом в работах Б. П. Рудого [1]. В этом подходе модель МС (более широко - узлового сечения стыка элементов трубопровода) строится на обобщении задачи о РПР, при решении которой простые волны любого знака рассматриваются в изоэнтропном приближении - по выражениям вида (8) и (9), а для узлового сечения всегда предполагается возможность применения в наборе «замыкающих» соотношений его универсальной статической характеристики, заданной одним или более критериальными уравнениями. Так, для излагаемых ниже моделей МС, относящихся к классу «скачков сечения», необходимым «замыкающим» соотношением может служить зависимость вида (1) или (при отсутствии данных о влиянии числа Рейнольдса на потери полного давления в «сжимаемом» течении) - вида (2).
Рассмотрим РПР на МС, установленном между участками каналов, имеющих в общем случае неодинаковое и переменное по длине сечение. Исходные данные задачи - параметры од-
номерного потока, свойств газовой смеси и площади граничных сечений примыкающих каналов рь Ть М1, Я1, уь ^ и р5, Т5, М5, Я5, у5 (рис. 1).
Рис. 1. Течение при распаде разрыве на скачке сечения (модель МС типа ДИАФРАГМА)
Из решения автомодельной по х/1 задачи о РПР могут быть найдены параметры потока в зонах 4 и 3, а по ним - искомые плотности потоков массы, импульса и энергии на границах крайних расчетных ячеек, необходимые для обновления решения в них на текущем временном шаге. Описываемая далее расчетная процедура использует приближение смеси совершенных газов; она используется в СИМ «Альбея» в качестве модели течения через МС типа ДИАФРАГМА.
В процедуре расчета РПР вначале определяется направление течения, для чего рассчитываются параметры нестационарного торможения потоков в зонах 5 и 1:
Р5 = Р5 /Р (М5 , У5 ), Р1 = Р1 /Р (М1, У1), ТГ= Т5/?(М5,У5), Т5 = Т1/(М1, У1),
= С5/а (М5, у5), с = с /а ' М, У1).
Этот (рекомендуемый в [1]) прием интерпретируется как появление в сечении МС на неопределенно малое время тонкой заслонки (рис. 1), после исчезновения которой, в силу обратимости торможения и разгона в простых изоэнтропных волнах, разовьется волновая структура течения с теми же значениями искомых параметров Р4, Т4, М4, Я4, у4, Р3, Т3, М3, Я3, у3, как в отсутствие заслонки.
Газ при РПР будет двигаться в направлении меньшего давления нестационарно заторможенного потока. Так, при Р5 >Р1 картина течения будет соответствовать показанной на рис. 1: течение в зонах 4, 3 и 2 вправо М4>0, М2>0, У3 = у4 = у5, Я3 = Я4 = Я5, 71 = 72, ^2 = - с образованием контактной поверхности (КП),
разделяющей газы разного «происхождения». В противном случае (р5 <Р1 ) направление течения будет обратным, но расчетная схема (рис. 1) и процедура останутся в силе после переиндексации параметров в исходных данных.
После выражения соотношений на всех элементах волновой картины течения (рис. 1) при РПР через газодинамические функции получим систему уравнений, связывающих неизвестные М4, М3, и М2:
1
р5 = Р(М2, ТО Р(М4, у5) _
Р1 р" 4, у 5 ) У5) а43 , ...)
< = М2 а' (М2, у1) а(М4, у5)
(16)
.(17)
с1 Мз а"(М4, у5) а(Мз, у5)" Ч(-М4, у5)-Р = 9(Мз, у5)^1а4з(М4, ...). (18)
В уравнение (18) входит ГДФ д(М у), связывающая расход О с параметрами стационарного торможения в том же сечении, см. напр., (19).
Решение задачи о РПР, описываемое системой (16)-(18), удобно находить, уточняя одну из неизвестных - например, М4, в интервале (0, М4тах], где М4тах определяется зависимостью1
ММ* ■■■).
Полученное на очередной итерации значение М4 определяет статические параметры и параметры стационарного торможения в зоне 4:
Р4 = Р5Р"(М4, у5), Г4 = Т"т"(М4, у5 ),
Р* = Р1/Р(М4, у5 ) , Т4 = 74/Т(М4, у5 ) .
Расход и параметры торможения в зоне 3 за МС определяются с использованием всех трех условий
О3 = О4 = т^(М4,^рР* ,
73* = 74*, Р* = Р4а43(М4).
Число М3 находится из уравнения сохранения массы и энергии в зоне 3
О,
= „, д(М з, у5)Р*
: т
(19)
после чего определяются статические парамет-
1 Единственность решения обеспечивается тем, что такая зависимость, полученная из аналитической модели МС, или численным расчетом течения на нем, или же продувкой МС, монотонно убывает; причем (что для обычных условий вполне корректно) можно считать М4 < 1. В дальнейшем, несколько упрощая, предполагаем, что зависимости для а задаются функциями одного числа М (здесь: а43(М4)), а решения задач о РПР располагаются на ветви, где а монотонно убывают. В противном случае решение найдется по модифицированной процедуре в точке (а< а(М), М = Мтах), т. е. на вертикальном участке характеристики потерь, и будет соответствовать режиму «запирания».
ры состояния в этой зоне и скорость потока:
7з = Т3т(М 3, у 5), Рз = Р*зР.(М 3, у 5), с3 =^]у5Я5Т3, и3 = М3с3.
Затем из соотношения Р2 = Р3 = р: п(М2, у^ определяется число М2, а по нему - остальные параметры потока в зоне 2:
С2 = с[а (М2, у!), и2 = М2С2.
Итерации по М4 прекращаются по достижении равенства и2 = и3 с заданной точностью.
Полученное решение удовлетворяет как соотношениям на фронтах простых волн и на КП, так и стационарным соотношениям, описывающим течение через МС как адиабатное и необратимое. Потоки массы О, и энергии (ОН )■ в примыкающих к МС зонах 4 и 3 одинаковы, потоки импульса вычисляются из выражения (Ои+рР), (■ = 3, 4 - индексы сечений или зон).
Если же с одной стороны от МС расположена емкость, а не канал, должна решаться задача о РПР иного вида, когда или
(рис. 2). Описываемые далее расчетные процедуры также используют приближение смеси совершенных газов и реализованы в СИМ «Аль-бея» в качестве моделей течения через МС типа КЛАПАН. При РПР на МС типа КЛАПАН возможны два направления течения - с истечением из емкости (рис. 2, а), которое реализуется при р0>р: , и с втеканием в емкость (рис. 2, б), реализуемом при Р0<Р\ . Если же Ро~Р\ (или сечение МС перекрыто запорным органом), принимается О3 = (ОН )3 = 0, а поток импульса в канал - (Ои+рР)3 = р, р («тривиальный» РПР на МС).
При РПР с истечением из емкости (р0>р: , рис. 2, а) КП отделяет поступающий в канал газ от газа в трубопроводе. Режим стационарного перетекания через МС задается отношениями Т3 = Т0 = Т0 иР3 = Р0о0(М3, ■..). Система уравнений модели данного МС на данном режиме течения связывает искомые числа М3 и М2 через ГДФ стационарного и нестационарного торможения:
1
Р0 = Р (М 2, у:)__
р! лМ у 0) а0зМ ...):
с^ = а (М2, у,) _М2 с] а(М3, у0) М3.
(20) (21)
Решение задачи о РПР по данной модели также удобно находить итерационным методом, уточняя, например, М3 в интервале (0, М3тах], определяемом зависимостью, которая задается
чаще всего как о-03(М3). Тогда для М3 в очередном приближении рассчитываются параметры в зоне 3
Р3 = Р0°03 (М3 )р(М3 , Уо X Т3 = Т)^(М3, Уо ), а также скорость потока
«3 = М3С3 = М3У1 УоЯТ3 .
Далее, аналогично тому, как описано выше для РПР на МС ДИАФРАГМА, определяются параметры в зоне 2; условие и3 = и2 также может служить для проверки сходимости итераций.
Рис. 2. Два режима течения при РПР а МС КЛАПАН: а — истечение из емкости; б — втекание в емкость
При РПР с втеканием в емкость (Р0<Р1, рис. 2, б) согласно принятой модели в канале формируется простая изоэнтропная волна и для расчета течения в зоне 3 нужно найти М3 из уравнения
= р \М 3, У1)
Р1 р(М3, У1)
• о
30
(М3|, ...).
В (22) чаще всего подставляется характеристика потерь на МС для режима втекания в емкость вида а30 = ст30(М3), здесь ^30 = Р0 /Р3 = Р0/Р3 . Поиск корня М3 уравнения (22) в интервале (0, |М3тах|] должен проводиться в общем случае также итерационно.
3. МОДЕЛЬ ПУТЕВЫХ ПОТЕРЬ
Источниковый член в уравнении количества движения (4) в нашей модели «замыкается» через коэффициент трения X, для вычисления локального значения которого имеется множество полу-эмпирических зависимостей; в случае газа они должны соответствовать общему критериальному уравнению
А = А(яе, А, М, у, Рг, в),
где А = А3э - относительная эквивалентная высота неровностей для шероховатой трубы и т. д. При умеренных дозвуковых значениях М в сечениях трубы не нужно учитывать сжимаемость, т. е. достаточно использовать зависимости вида А = А(яе, а) . Во всех расчетах данной работы было принято Д=0 и для развитого турбулентного режима течения в трубе применялась формула, рекомендованная в [9] для диапа-
зона Ф103<Яе<Ф106:
= 2^10 (яеЛ)-0,8.
л/А
(23)
Локальная плотность потока импульса на стенке трубы в (4) через X выражается как
ри|и| ^ ТУ = -А—— —. " 23э П
Локальная плотность потока энергии от теплоотдачи в стенку канала может быть оценена при помощи аналогии Рейнольдса, тогда указанная величина, фигурирующая в уравнении энергии (5), определяется выражением
Чу
2Л
Т - Т * £
' ™ П
(25)
Формула (25) получена в предположении о том, что «эффективное» число Ргэфф - во всем сечении равно 1. Поскольку как «молекулярные» числа Рг, так и «турбулентные» Ргт в потоке газа не достигают единицы, оценка теплового потока по (25) будет завышенной. В нашей модели в (25) введен поправочный коэффициент Рг0,5, рекомендуемый в одной из моделей теплоотдачи для турбулентного пограничного слоя; для двухатомных газов и их смесей (Рг^0,72) значение этого коэффициента 0,72°,5~0,85.
4. МЕТОДИКА И РЕЗУЛЬТАТЫ ПРОДУВОК МЕСТНЫХ СОПРОТИВЛЕНИЙ
Необходимые для расчетов движения волн по трубопроводу установки характеристики потерь полного давления на его местных сопротивлениях (МС) определены статическими продувками на специальном стенде (рис. 3).
с
Рис. 3. Продувочный стенд
В качестве МС в экспериментах использованы круглые шайбы с номинальными диаметрами отверстий = 18 мм, 14 мм, 10 мм и 6 мм, устанавливаемые как между двумя участками трубопровода, так и на его конце. Для фиксации шайб на конец трубы наворачивалась специальная гайка с отверстием = 24 мм, при установке которой без шайбы дополнительно получалось МС с указанным номинальным диаметром. Для получения наименьшего сопротивления потоку при втекании в трубу была изготовлена «гайка» специальной формы, часть поверхности которой была выполнена по лемнискате Бер-нулли (ЛБ) и чисто обработана (рис. 4).
Рис. 4. Отрезки труб, шайбы, лемниската Бернулли и фиксирующая гайка
Потери полного давления определялись обработкой измерении 003, 030 043 для различных значений числа Мт для всех вариантов установки сопротивлений на трубопровод по приведенной ниже методике. Для случая «свободного» втекания в емкость из трубы (через гайку без шайбы: = 24 мм и через ЛБ) характеристика потерь может быть задана как а30 = п(Мт, у = 1,4), что эквивалентно условию Рт = Р0
на срезе трубы для Мт < 1,0.
Измерения проведены в потоке атмосферного воздуха (рис. 5), расход которого G создавался насосом 1, а измерялся расходомерами 2 объемного типа РГ-250 и РГ-40 (последний применялся при величине объемного расхода менее 40 м3/час). Давление p2 в ресивере-успокоителе 3 измерялось образцовым вакуумметром 4 класса точности 0,25. Для разрежений Ap02 = p0-p2 < < 1,2 м вод. ст. применялся водяной пьезометр, не показанный на рис. 5, им же контролировался перепад давлений Ap01 на расходомере.
Атмосферное давление p0 в лаборатории определялось по настенному барометру-анероиду (±0,5 мм рт. ст.), температура T0 - по лабораторному термометру (±0,25 °C). Интервалы времени (20 - 75 с) для измерения G отмечались по секундомеру «Агат» (±0,2 с). Измерение для каждой экспериментальной точки производилось трижды.
При обработке результатов измерений в методике использовано допущение об адиабатно-сти течения T = T0. Параметры потока приводились к сечениям, примыкающим к МС, для чего выполнялся формальный учет трения о стенку на длине труб (l1 = 1038 мм, l2 = 931 мм), которые считались цилиндрическими с dx = 24,1 мм.
Так, для измерений по схеме на рис. 5, а число M2 в выходном сечении трубы определялось из уравнения
G = my(M2'Д^p2,
VT7
где m = 0,0404 м'1-с-К0,5 для воздуха с у = 1,4 и R = 287,1 Дж/(кгК), y(M, у) = q(M у)/я(М у), см. (15) и, например, (19). По известным значениям G, M2, p2 с учетом T = consti и G = const2 численным решением ОДУ
dJ dx
полученного из (11) при dF/dx=0 и П = ndx для
' *
участка длиной l1 вычисляются px и Мт, приведенные к сечению за МС, после чего а03 (для данного режима течения по MT) берется по определяющему соотношению px/p1. Таким способом определена и зависимость а03ЛБ (Мт) для ЛБ, использованная для обработки измерений по схемам на рис. 5, б и в.
Измерения по схеме на рис. 5, б обрабатывались несколько иначе. Решением уравнения
G = m
q(M, g)FT pxs
03 ЛБ (M)
л/T*
определялись параметры М, р = р1с03ЛБ(М) на входе в трубопровод, а после расчета течения с трением по (24) и (23) - параметры рт и Мт перед шайбой, установленной на выходе из трубы, после чего а03 определялся как рт /р2.
Рис. 5. Схема воздушного продувочного стенда для снятия характеристик: а - аоз(Мт); б - озо(Мт); в - о^зМ
1
0,9
8 0,8 ö
0,7
0,1 0,2 Мт 0,3 0,4 0,5
1 \ 1 \ 1
0,9 \ 0,9 \ 0,9
1- 0,8 Л 0,8 \ 0,8
со <8 0,7 ■ \ 0,7 \ 0,7
0,6 0,5 ж Г 0,6 0,5 х Д ! 0,6 0,5
0 0,3 Мт
0 0,16 Мт
0 0,06 Мт
Рис. 6. Зависимости с03(Мт) на режиме втекания в трубопровод: а - ё0 = 24 мм (ЛБ); б - ё0 = 24 мм; в -ё0 = 18 мм; г - ё0 = 14 мм; д - ё0 = 10 мм; е -ё0 = 6 мм
Процедура обработки измерений по схеме на рис. 5, в представляет собой комбинацию двух описанных. В этом случае определяемая для каждого измерения величина о43 вычисля-
лась как отношение полных давлений, приведенных к сечению шайбы, «режимным» параметром являлось число Мт перед шайбой. Полученные таблицы экспериментальных точек были аппроксимированы (рис. 6 - 8) зависимостями, которые при Мт ^0 стремятся к квадратичным.
1
0,9
0,8
§ 0,7
0,6 0,5
\ \ + а
\ \ хб \
\ \ О в \
\ а г \
0
0,1
Мт
0,2
0,3
Рис. 7. Зависимости ^30(Мт) на режиме истечения из трубопровода: а — d0 = 18 мм; б — d0 = 14 мм; в — d0 = 10 мм; г — d0 = 6 мм
0
0,1 Мт 0,2
0,3
Рис. 8. Зависимости ^43(Мт) для шайб в трубопроводе: а —d0 = 18 мм; б — d0 = 14 мм; в — d0 = 10 мм; г — d0 = 6 мм
5. ОДНОЦИКЛОВАЯ УСТАНОВКА И УСЛОВИЯ ЭКСПЕРИМЕНТОВ
Нестационарное движение ВКА по длинному неразветвленному трубопроводу с МС воспроизводилось в экспериментах на одноцикло-вой установке, содержащей генератор волн и присоединенный к нему трубопровод. Схема генератора волн представлена на рис. 9.
Основной элемент генератора волн - емкость 1, соединенная с трубопроводом клапаном 2, открываемым при помощи рессоры 3; подъем и время открытия клапана регулируются гайкой 4. Рессора удерживается во взведенном состоянии защелкой 5, соединенной с сердечником соленоида 6. При подаче тока на соленоид защелка освобождает рессору, которая приводит в движение клапан. После того, как зажимы 7 на
0
рессоре выходят из втулки 8 и освобождают клапан, он под действием возвратной пружины 9 закрывается.
Рис. 9. Схема генератора волн одноцикловой установки
К емкости присоединен трубопровод 10, по которому распространяются ВКА. Разрежение или повышенное давление в емкости создается либо компрессором, либо вакуумным насосом 12. Величина давления регулируется вентилем 13. Рессора вводится в зацепление с клапаном и защелкой при помощи рычага 14.
В измерительную систему входят: один или два датчика давления, усилитель-преобразователь, аналого-цифровой преобразователь (АЦП) и персональная ЭВМ. Первичный преобразователь измерительной системы - малоинерционный датчик давления с чувствительным элементом в виде кремниевой пластинки с нанесенным на нее активным сопротивлением, включенным по мостовой схеме. Вторичным преобразователем служил индуктивный высокочастотный многоканальный преобразователь (усилитель) 6705+ фирмы AVL. Измеряемый сигнал с выходов преобразователя подавался на АЦП PCI-1711 фирмы Advantech, установленный в слот PCI материнской платы ПЭВМ. Данный АЦП с разрядностью 12 бит обеспечивал преобразование аналогового сигнала в диапазоне напряжений от -10 В до +10 В в цифровой эквивалент (отсчет) от 0 до 4095. Взаимодействие с АЦП обеспечивалось специально разработанной программой для ЭВМ, позволяющей выполнять тарировки измерительных каналов, воспроизводить записи сигналов датчиков на экране и записывать полученные данные в файл.
Перед проведением серии измерений выполнена тарировка обоих измерительных каналов с датчиками давления. При тарировке на датчики подавались разрежения от 0 до -0,64 кГс/см c шагом в 0,16 кГс/см2 (и такие же
уровни избыточных давлений) с применением образцовых манометра и вакуумметра класса точности 0,4 и шкалами на 1 кГс/см2. Тариро-вочные зависимости каналов для учета их нелинейности аппроксимируются программой обработки по методу наименьших квадратов многочленами третьей степени:
Уу = рХу) = Ах3 + Вгх) + Сгх} + Ц, I = 1, 2,
где Ху - величина отсчета АЦП.
Перед тарировкой и перед каждым измерением показания датчика (датчиков) при поданном на них атмосферном давлении приводились к нулю - для исключения его дрейфа. Дискретность по времени считывания показания в измерительном канале составляла 40 мкс. Приведенные ниже экспериментальные кривые подвергнуты осреднению по трем точкам с целью фильтрации помех: <у> = (уу+1+ уу+ уу-1)/3.
Трубопровод экспериментальной установки (рис. 10) собирался из труб (материал - нержавеющая сталь), в которых с высокой точностью выдержана круглая цилиндрическая форма внутренней поверхности с шероховатостью не более ^10.
1
Рис. 10. Вид одноцикловой установки с трубопроводом и измерительной системой
В варианте исполнения, показанном на рис. 11, а - лемниската Бернулли (ЛБ) или шайба устанавливалась на конце трубопровода, в варианте на рис. 11, б шайбы с различным номинальным диаметром устанавливались в месте соединения двух труб примерно одинаковой длины.
Таким образом, в проведенном расчетно-экспериментальном исследовании прохождения ВКА через шайбы вариант а может служить для выявления погрешностей модели МС КЛАПАН, вариант б - модели МС ДИАФРАГМА.
В качестве элементов трубопровода, примыкающих к шайбам или ЛБ, использованы отрезки труб, примененные при статических продувках этих МС.
Суммарные длины участков трубопровода: 11 = 3028 мм, 12 = 3692 мм, диаметры -31~24,1 мм, ё2 = 24,1 - 24,8 мм (реальные профили составных труб учтены в расчетах), координаты датчиков: /д1 = 781 мм и 1д2 = 992 мм (см. рис. 11).
Рис. 11. Схема двух вариантов исполнения трубопровода: а — с МС на конце трубы; б — с МС на стыке труб
Объем емкости генератора волн определен расчетным путем на основе поэлементного обмера: V = 748,2 ± 5,7 см3. Исходными ВКА во всех экспериментах были волны разрежения.
6. СРАВНЕНИЕ РАССЧИТАННЫХ И ИЗМЕРЕННЫХ р(0
Для проверки изложенных выше моделей и методик, с их применением выполнены численные расчеты нестационарных течений в трубопроводе для условий выполненных экспериментов. Для этого в системе имитационного моделирования (СИМ) «Альбея» (рис. 12) были построены модели обоих вариантов исполнения трубопровода (рис. 11), где задавались их действительные размеры, свойства и начальные параметры рабочего тела (атмосферного воздуха), а также полученные из статических продувок зависимости для потерь полного давления на шайбах и лемнискате Бернулли (рис. 6 - 8).
Течение через клапан генератора волн одно-цикловой установки рассчитывалось в СИМ «Альбея» также по модели РПР на МС.
В исходных данных задавались зависимости для коэффициента восстановления полного давления а03 и а30 от Мт в примыкающем сечении трубы, полученные расчетами для нескольких значений высоты его подъема к^ по приближенной методике, описанной в [2, 10]. На рис. 13 показаны графики полученных зависимостей а30(Мт, к) для случая течения газа с у = 1,4 в направлении из канала в емкость.
1
0,8 0,6
о 0,4
со
о
0,2 0
0
0,2 Мт 0,4
0,6
Рис. 13. Графики расчетных зависимостей а30(МТ) для клапана одноцикловой установки для разных значений высоты подъема клапана: к = 1,5;
3,0; 4,5; 6,0; 7,5 мм
В расчетах волнового течения в трубопроводах применено уравнение состояния для воздуха как смеси идеальных газов с теплоемко-стями, зависящими от Т. Соотношениями (24) и (23) учтено локальное напряжение от трения на стенке; ввиду того, что при движении ВКА по трубе Т = (х, /) Ф Ту = Т0, в расчетах учтен локальный тепловой поток от стенки - соотношением (25) с поправочным коэффициентом.
Применение для расчета характеристики клапана, полученной по приближенной методике, оправдано тем, что амплитуда и форма исходной ВКА разрежения в расчетах задавалась путем подбора закона изменения единственной зависимости, выступающей в роли граничного условия в сечении трубопровода, соответствующем седлу клапана. Однако поскольку в СИМ «Альбея» возможность задания граничных условий, например, вида р(/), не реализована, подбором определялся закон открытия клапана к(/) генератора волн.
Зависимость к(/) подбиралась с помощью программной утилиты многопараметрической оптимизации, входящей в СИМ «Альбея».
Рис. 12. Модель трубопровода в СИМ «Альбея»
Рис. 14. Давление на Дь = 24 мм (ЛБ)
Для каждого расчета определялись наборы оптимальных значений 16 параметров закона к(/) открытия клапана с учетом двух подскоков после его первой посадки. Минимизировано среднеквадратическое отклонение расчетного давления от измеренного на участке исходной волны. Дальнейшее расхождение расчетной и измеренной кривых давления должны объясняться погрешностями расчета по модели, которые вызваны схематизацией в описании процессов как на участках между МС, так и при прохождении самих МС.
Рис. 15. Давление на Дь = 24 мм (шайба)
Расчеты проведены с шагом по времени & = 10 мкс, применен двухэтапный (имеющий второй порядок аппроксимации по /) явный метод обновления искомых переменных в ячейках расчетной сетки.
Рис. 16. Давление на Дь d0 = 18 мм
Рис. 17. Давление на Д1; d0 = 14 мм
В модели ТРУБКА использована дискретизация повышенной точности [5] по пространственной переменной x, что (при 460 и 560 ячейках соответственно в 1-й и во 2-й трубках) гарантирует более чем достаточную для практики сходимость численного решения задачи к точному. Расчет процесса в трубопроводе на рис. 11, б для интервала tj = 0,2 с модельного времени занимал 82 секунды процессорного времени одного ядра Intel Core 2 Duo при тактовой частоте 1,5 ГГц.
На рис. 14 - 19 приведены расчетные и измеренные зависимости от времени для давления в сечении Д1 трубопровода вида рис. 11, а при установленных на свободном конце трубопровода соответственно ЛБ (рис. 14) и шайбах с номинальными диаметрами d0 = 24, 18, 14, 10 и 6 мм. Ход измеренного давления показан прерывистой линией.
Рис. 18. Давление на Д1; й0 = 10 мм
с
Рис. 19. Давление на Д1; й00 = 6 мм
Рис. 20. Давление на Д1 и Д2 й0 = 24 мм (шайбы нет)
Рис. 22. Давление на Д1 и Д2 й0 = 14 мм
Рис. 23. Давление на Д1 и Д2 й0 = 10 мм
На рис. 14 видно, что избыточное давление в волне сжатия, отраженной от ЛБ (минимальные потери давления при втекании), превышает перепад давления в исходной ВКА разрежения. Заметно снижение интенсивности отраженной волны при переходе от ЛБ к шайбе с й0 = 24 мм (рис. 15). Напротив, при й0 = 6 мм отражение от МС на конце трубы происходит почти как от закрытого конца - отраженная ВКА имеет тот же знак (разрежение) и примерно такой уровень минимального давления, что и исходная (рис. 19).
Рис. 24. Давление на Д1 и Д2 й0 = 6 мм
Рис. 21. Давление на Д1 и Д2 ё0 = 18 мм
Во всех расчетах этой серии отмечается близкое совпадение кривых давления в части формы и амплитуды, по крайней мере, на участках первых отраженных волн; на расчетных кривых давления на протяжении всего интервала модельного времени отсутствуют какие-либо нефизичные дефекты.
Следует отметить характерное для всех экспериментов серии нарастание погрешности расчета протекания кривых давления на датчиках по времени, а именно - занижение расчетом характерного периода волнового процесса. Погрешность времени регистрации характерных участков волн датчиками нарастает приблизительно линейно, и это нарастание, по-видимому, слабо связано как с интенсивностью так и со знаком волн. Такого порядка временное рассогласование (-2,0 ... -2,8 %, по рис. 14) моментов прибытия волн к датчикам отмечено автором еще в [2] на основе данных, полученных с другой измерительной системой и обработанных другой расчетной программой; это практически исключает вероятность «промаха» от задержек в измерительной системе, ошибочного задания длин труб или уравнения состояния, искажающего величину скорости звука. Обнаруживаемое расхождение, скорее всего, следует из неадекватности квазиодномерной модели применяемого вида для нестационарных процессов даже в случае из протяженных (¡/ё — 140) участков.
На рис. 20 - 24 сопоставлены зависимости рг(0 в сечениях Д1 и Д2 для трубопровода вида рис. 11, б с номинальными ё0 = 24, 18, 14, 10 и 6 мм соответственно; во всех случаях на свободном конце трубопровода (суммарной длиной ¡1 +¡2—6,7 м) устанавливалась ЛБ.
В данной серии также не наблюдается каких-либо нефизичных дефектов на расчетных кривых давления («всплески» на экспериментальных кривых - рис. 14 - 16 и рис. 20 - 22 -помехи в электрической цепи соленоида в момент отпускании кнопки ПУСК). Согласование величин расчетного и измеренного давлений в данной серии несколько хуже - что видно даже по р(0 на Д2 в отсутствие МС на стыке каналов (рис. 20); это, возможно, проявление погрешности в учете путевых потерь при большей длине трубопровода.
Погрешность, проявляющаяся в уменьшении периода волнового процесса, характерна и для расчетов данной серии (рис. 20 - 24).
ЗАКЛЮЧЕНИЕ
Расчетно-экспериментальное исследование, предпринятое для проверки расчетной методи-
ки, построенной на одномерной модели нестационарного течения применительно к протяженным неразветвленным трубопроводам, подтвердило в целом высокую достоверность результатов расчета.
Действительно, оценивая достоверность результатов расчета величиной относительной погрешности определения амплитудного значения давления в первой отраженной от МС волне (а именно, относя абсолютную погрешность к возможному разбросу амплитудного давления для отраженной волны любого знака) получим, что для рис. 14 - 24 эта величина не превышает 1,6%.
Такой уровень погрешности является значимым при принятом методе измерений; он может объясняться совместным влиянием погрешностей задания амплитуды исходной волны и определения характеристик потерь на МС, а также погрешностями, вызванными принятым подходом к учету путевых потерь.
В дальнейшем, со снижением амплитуд волн, наблюдается нарастание абсолютной (и тем более - относительной) погрешности в определении величины давления расчетом.
Следует отметить, что полученные численные решения свободны от нефизичных дефектов.
Поэтому можно (с оговоркой, см. ниже) утверждать, что в данном исследовании показана достаточно высокая степень достоверности результатов моделирования нестационарных течений газа в протяженных неразветвленных трубопроводах по методике, включающей
• использование в качестве модели системы уравнений, описывающих течение в канале в квазиодномерном приближении;
• «замыкание» данной модели эмпирическими соотношениями для локальной интенсивности трения и теплоотдачи в стенку канала, справедливыми для стационарных условияй;
• использование моделей нестационарного течения газа через МС, основанных на обобщении задачи о распаде разрыва и замыкаемых зависимостями вида о(Мт) для описания течения непосредственно на МС;
• методику проведения и обработки статических продувок МС для получения указанных характеристик МС;
• применение явного консервативного метода для решения уравнений газовой динамики.
Подтвержден обнаруженный ранее [2] эффект завышения расчетом в квазиодномерном приближении скорости распространения ВКА по трубопроводу, снижающий точность моделирования. Наиболее вероятная причина эффекта -
неадекватность системы уравнений «квазиодномерного» приближения примененного вида для описания волнового движения в канале.
Можно предположить, что обусловленная учетом неоднородности распределения параметров потока в сечении канала коррекция уравнений, приводящая к изменению выражений для наклона характеристических кривых (7), способна в значительной степени компенсировать данный эффект.
СПИСОК ЛИТЕРАТУРЫ
1. Рудой, Б. П. Прикладная нестационарная гидрогазодинамика : учебн. пособие / Б. П. Рудой. Уфа : УАИ, 1988. 184 с.
2. Черноусов, А. А. Определение гидравлических характеристик местных сопротивлений в газовоздушных трактах ДВС вычислительным экспериментом: дисс... канд. техн. наук /
A. А. Черноусов. Уфа : УГАТУ, 1998. 164 с.
3. Черноусов, А. А. Консервативная сеточно-характеристическая схема для расчета нестационарных течений в трубопроводах / А. А. Черноусов // Вопр. теории и расч. рабочих проц. тепл. дв-лей. Вып. 21. Уфа : УГАТУ, 2007. С. 246-254.
4. Горбачев, В. Г. Система имитационного моделирования «Альбея» (ядро). Руководство пользователя. Руководство программиста : Учеб. пособие /
B. Г. Горбачев, С. А. Загайко, Н. В. Рудая и др. Уфа : УГАТУ, 1995. 112 с.
5. Черноусов, А. А. К выбору метода численного интегрирования уравнений одномерного движения газов в каналах ДВС / А. А. Черноусов // Известия вузов. Машиностроение. 2007. № 6. С. 34-46.
6. Рождественский, Б. Л. Системы квазилинейных уравнений и их приложения к газовой динамике / Б. Л. Рождественский, Н. Н. Яненко. М. : Наука, 1978. 687 с.
7. Самарский, А. А. Разностные методы решения задач газовой динамики / А. А. Самарский, Ю. П. Попов. М. : Наука, 1980. 352 с.
8. Дулов, В. Г. Распад произвольного разрыва параметров газа на скачке площади сечения / В. Г. Дулов // Вестник ЛГУ. Сер. матем., мех. и астрон. 1958. № 19. Вып. 4. С. 76-99.
9. Абрамович, Г. Н. Прикладная газовая динамика: учебн. руководство / Г. Н. Абрамович. Ч. 1. М. : Наука, 1991. 600 с.
10. Рудой, Б. П. Определение коэффициентов гидравлических потерь в вычислительном эксперименте / Б. П. Рудой, А. А. Черноусов // Актуальные проблемы авиадвигателестроения. Уфа : УГАТУ, 1998. С. 189-197.
ОБ АВТОРЕ
Черноусов Андрей Александрович, доц. каф. ДВС УГАТУ. Дипл. инж.-мех. по ДВС (1994 г., УГАТУ). Канд. техн. наук по тепловым двигателям (1998 г., УГАТУ). Иссл. в обл. числ. моде-лир. газообмена и раб. проц. ДВС, выч. гидрогазодинамики.
