CC BY
42
ФЕНОМЕНОЛОГИЧЕСКАЯ МОДЕЛЬ НЕСТАЦИОНАРНОГО ТЕЧЕНИЯ СЛАБОСЖИМАЕМОЙ ЖИДКОСТИ ЧЕРЕЗ ДРОССЕЛЬ
В. Б. Синильщиков
Балтийский государственный университет «Военмех» им. Д. Ф. Устинова, доцент, [email protected]
1. Общая характеристика нестационарных задач гидравлики и методов их решения. Задача о нестационарном течении сжимаемой жидкости в трубопроводах и каналах, в том числе имеющих местные сопротивления, такие как сужения, расширения, дроссели, клапаны и др., имеет важнейшее практическое значение. Можно отметить такие ее приложения, как работа гидродемпферов при ударном и высокочастотном нагружении амортизированных систем, нестационарные (в т. ч. автоколебательные) режимы работы клапанов, гидроудар в сложных и разветвленных гидросистемах и некоторые другие. Универсальные методики расчета подобных процессов до настоящего времени отсутствуют, и на практике (например, при выборе жесткости клапана и сечения подводящего канала, параметров и расположения компенсаторов и гасителей гидроудара) руководствуются, в основном, опытом проектирования или имеющимися экспериментальными данными. Основной проблемой, затрудняющей решение подобных задач, является необходимость одновременного и совместного учета двух групп процессов: гидродинамических (связанных с установившимся или близким к нему течением жидкости) и волновых. Первая группа процессов достаточно хорошо описывается обобщенным уравнением Бернулли для потерь полного напора между характерными сечениями канала (I и II), имеющим вид (см. [1])
Здесь n и m — число местных сопротивлений и участков постоянного сечения между сечениями I и II, щ — средние скорости в сечениях (из условия несжимаемости: UjFj = ukFk = const), pi —давления, hi —высота, Fk, Щ и lk —площади и периметры поперечных сечений и длины участков, Tk — среднее значение касательных напряжений на стенке k-го участка канала, Qi — коэффициенты местных потерь, р — плотность жидкости. Коэффициенты ai учитывают различие между средней скоростью, входящей в выражение для расхода, и среднеквадратической, которая определяет скоростной напор. Отметим, что эти коэффициенты достаточно близки к единице, причем стандартные рекомендации по их определению относятся к установившимся режимам течения и не могут быть использованы для нестационарного случая. Поэтому при записи уравнений, описывающих нестационарные режимы течения, обычно принимают i = 1. Достаточно полный набор экспериментальных данных и аппроксимаций величин Q и Tk для стационарного случая приведен в [2].
© B. Б. Синильщиков, 2011
Другой группой процессов, которые необходимо рассматривать при решении подобных задач, являются волновые (акустические). Они связаны с распространением по каналам (трубам) волн сжатия-разрежения. Подобные процессы также достаточно хорошо изучены. Для канала постоянного сечения они описываются уравнениями, полученными Н. Е. Жуковским [3]:
ди
рт =
др
дх
тП
~У'
др о ди
т = ~рс~ді-
(2)
Здесь с — скорость звука (с учетом податливости стенок канала и наличия нераство-ренного воздуха [1]). Система (2) отличается от обычных уравнений одномерной акустики учетом трения на стенке. На основе данной модели в [4] рассматривается целый ряд задач о гидроударе в трубах. Для их решения используются аналитические методы (метод характеристик, метод Фурье). В современных работах для решения подобных задач используются численные методы: метод Годунова, метод конечных элементов [5]. Вместе с тем остается до конца нерешенной проблема стыковки результатов, получаемых на основе соотношений (1) и уравнений (2). Для участков, имеющих постоянное сечение и не содержащих местные сопротивления, система (2) в квазистационарном случае эквивалентна (1) при а = 1. Отдельную проблему представляет определение напряжений трения при нестационарных режимах. Подробные исследования этой проблемы проведены в работах Д. Н. Попова [6-7]. Как показано в этих работах, нестационарное значение коэффициента сопротивления на стенке может существенно отличаться от стационарного. Ниже рассмотрим проблему нестационарных режимов течения через местные сопротивления на примере дросселя.
2. Задача нестационарного течения жидкости через дроссель. Существующие модели. Дроссель — рис. 1—один из типовых элементов гидросистем. Здесь и ниже индексами «1» и «2» обозначены параметры в зонах одномерного течения в каналах (камерах) перед дросселем и после него, а индексом «12»—в самом дросселе. Для определенности рассмотрим случай, когда торцевые стенки обоих каналов перпендикулярны осям каналов. Отметим, что такие местные сопротивления, как внезапное сужение или расширение и тонкое отверстие, могут рассматриваться как частные случаи дросселя при = ^2, и Ню = 0 соответственно.
Рис. 1. Схема течения жидкости через дроссель.
Сначала рассмотрим случай стационарного течения. С учетом малости размеров дросселя и допущения а = 1, уравнение (1) примет вид
Р1 + Р-^- ~ (р2 +Р~^) = (3)
В [2] приводятся выражения для коэффициента потерь напора С12 для разных режимов течения и сочетаний геометрических параметров. Обобщенную зависимость при
и 12 > 0 можно представить в виде
С12 = + (сч 1 - Л) + - /2) + (1 - /2)2) ё0йе + А7712. (4)
Здесь / = Дю/Д®, П12 = Л.12П12Д4Д12), коэффициент характеризует скругление кромок на входе (0 < ^' < 0, 5), к = к(п12) —влияние длины канала 1-2 на характер течения в канале (0 < к < 1, 35), ЛДв^) —сопротивление трения в канале 1-2, ^ = С^(Дв12, /1) и £^е = е^е(Дв12) — влияние вязкого сопротивления на входе и в канале при Дв12 < 105, Дв12 = 4^12и12/(П12^) —число Рейнольдса, V — кинематическая вязкость. В случае острых кромок и Дв12 > 105 при П12 > 2, 5 имеем £12 = 0, 5(1 — Д) + (1 — /2)2 + Л7712, а при 7712 -> 0 получим С12 = 0, 5(1 - Д) + 1, 35^/1 - /1(1 - /2)2 + (1 - /г)2-
Далее рассмотрим подходы, применяемые при решении нестационарных задач. Наиболее очевидным является подход, основанный на прямом решении осесимметричных или трехмерных уравнений гидродинамики сжимаемой жидкости численными методами. Однако опыт расчетов, проведенных автором с использованием метода Годунова, показывает, что при малой сжимаемости жидкости влияние схемной вязкости перекоса приводит к завышению расчетного значения сопротивления дросселя по сравнению с табличными значениями [2]. Для достижения приемлемой точности приходится значительно уменьшать размер ячеек (к примеру, чтобы снизить погрешность определения средней скорости истечения через дроссель при заданном перепаде до 5-10%, число ячеек на радиус минимального сечения дросселя должно составлять несколько десятков. В результате такая программа будет иметь крайне низкое быстродействие, что затрудняет ее использование для практических расчетов сложных систем.
Более эффективными являются подходы, основанные на одномерных моделях с использованием дополнительных сопрягающих условий на дросселе. Так, в [4] исследуется распространение волн сжатия в сложном трубопроводе, состоящем из участков, имеющих разные сечения. Границами участков являются внезапные сужения и расширения. На них ставятся сопрягающие граничные условия равенства расхода и давления. Таким образом, при подобном «акустическом» подходе давление на входе в сужающуюся часть, на торцевой стенке, а также в зонах одномерного течения в широком и узком каналах вблизи сужения считается одинаковым, т. е. реальное распределение давления не учитывается. Можно показать, что данное условие эквивалентно соотношению С12 = (Р/Д1)2 — 1, т. е. коэффициент потерь при сужении для установившегося течения получается отрицательным.
В большинстве более поздних работ [1, 5, 6] вместо условий равенства давления для сопряжения систем (2) на границах участков и иных местных сопротивлениях используются стационарные уравнения (3). Другим условием по-прежнему является условие равенства расхода в сопрягаемых каналах. В [6] отмечается, что «неустановив-шиеся течения в местных сопротивлениях еще мало изучены» и «вследствие этого при расчетах используют квазистационарные коэффициенты местных сопротивлений». В этом случае при установлении течения решение задачи на основе системы уравнений (2) и описанных условий сопряжения автоматически стремится к стационарному соотношению (1) при а = 1. Вместе с тем, такая модель не описывает нестационарную стадию.
В качестве примера рассмотрим задачу о распространении волны сжатия по длинной жесткой трубе с площадью сечения р = р2 = Д (/1 = /2 = /), в которой установлена дроссельная шайба с площадью сечения Д12 — рис. 2. Примем, что до прохождения волны жидкость покоилась.
Рис. 2. Схема взаимодействия волны сжатия с дроссельной шайбой: а-дии взаимодействия; г — линии тока при установившемся течении.
последовательные ста-
При скорости жидкости за волной ио перепад давления на фронте Дро = ри^е — рис. 2, а (здесь и ниже через Др обозначено избыточное по отношению к начальному давление). Часть фронта волны — рис. 2, б проходит сквозь дроссель, причем первоначальное значение скорости в дросселе равно ио. Другая часть фронта отражается от торцевой стенки шайбы. При этом давление в зоне отражения повышается и в полости 1 в районе входа в отверстие (на границе зон истечения и торможения) формируется радиальное течение, направленное к оси. Волна, прошедшая через зазор, расширяется до сечения канала за дросселем (давление при этом падает) и продолжает распространяться вправо, постепенно превращаясь в плоскую — рис. 2, в. Далее происходит постепенное формирование струи, истекающей из полости 1 в полость 2 — рис. 2, г. При этом скорость истечения постепенно возрастает, а перепад давления падает. Если труба настолько длинная, что прошедшая и отраженная волны, отразившись от ее концов, не успевают вернуться к дросселю до установления течения, то при установившемся течении скорость истечения, скорость и давление в зонах 1 и 2 определяются из уравнения (3), условия постоянства расхода и акустических соотношений Др1уст = ре(2ио — и1уст), Др2уст = реи2уст. Окончательные выражения для средних скоростей примут вид
2и0 и1уст /г\
^1уст — ^2уст — ■?» - - '
1 +
Если же, как в [5, 6], принять, что соотношение (3) справедливо с момента прохождения волны через дроссель, то и выражения (5) справедливы с этого момента. Это фактически означает, что волны мгновенно трансформируются в струйное течение, что явно противоречит здравому смыслу.
Задача определения времени установления течения через дроссель решалась в работе [7] в рамках модели несжимаемой жидкости. Для участка, ограниченного сечениями 1 и 2 (см. рис. 1), записывалось обобщенное уравнение Бернулли вида (1). Были получены выражения для суммарного коэффициента сопротивления ^12, подобные (4) в случае течения с образованием и без образования циркуляционной зоны в дросселе. Автор подробно проанализировал влияние нестационарных эффектов на коэффициент трения Л в канале 1-2 (в частности, была получена аналитическая зависимость от частоты, для гармонических колебаний малой амплитуды). Для других составляющих коэффициента сопротивления £12 использовались квазистационарные соотношения, хотя отмечалось, что нестационарные значения этих параметров могут отличаться от них. Наибольшую сложность представляет определение инерционного коэффициента, характеризующего приведенную массу струи, истекающей через отверстие. Для этого использовалась оценка, основанная на решении Кармана задачи о струйном обтекании
1 +
М0С12
с/2
12
I
пластинки ускоренным потоком жидкости, что, как показывает сравнение с результатами численных экспериментов, приводит к существенным погрешностям. Главный же недостаток такого подхода связан с неучетом волновой составляющей: при использовании такой модели расход через дроссель изменяется непрерывно. В частности, в рассматриваемой выше задаче (рис. 2) расход через характерное сечение дросселя после прохождения волны в соответствии с данной моделью равен нулю, а не иоР_2, как в реальности. В случае, если площадь сечения дросселя соизмерима с площадью канала, это различие весьма существенно.
3. Феноменологическая модель. Определение средних скоростей истечения. Как показано выше, в общем случае и волновые процессы и нестационарное течение после прохождения волн важны и должны учитываться при определении нестационарного расхода жидкости через дроссель. Вместе с тем анализ полей скоростей за фронтом волны (см. рис. 2, а-в) и при установившемся течении (см. рис. 2, г) показывает, что эти течения качественно отличны друг от друга (в частности, первое безотрывно, а второе связано с отрывом потока и образованием струи). Поэтому модели, в которых один из этих процессов не учитывается или рассматривается как частный случай другого, по мнению автора, не могут адекватно описать данную задачу. Так, не решает возникшие проблемы модификация рассмотренного выше подхода [7], основанного на обобщенном уравнении Бернулли, в соответствии с которой скорость за фронтом волны ио принимается в качестве начального условия для скорости истечения, а сам расчет начинается с момента прохождения волной сечения дросселя. Разные распределения скоростей за фронтом волны и при установившемся истечении приводят к разным инерционным коэффициентам при выражении кинетической энергии через скорость истечения, т. е. начальное значение кинетической энергии при таком подходе оказывается неверным. Другим недостатком такого подхода является сложность его обобщения на случай, когда акустические возмущения приходят к дросселю непрерывно.
Предлагаемый автором подход к решению подобных задач в какой-то мере является обобщением двух описанных выше. Он основан на разделении гидродинамических и акустических (волновых) процессов. Средняя скорость течения жидкости в характерном сечении канала 1-2 представляется в виде суммы двух составляющих:
и12 = и12ак + и12гндр* (6)
Первая из них (ниже будем называть ее акустической) характеризует составляющую средней скорости, которая изменяется мгновенно и связана с прохождением через дроссель акустических волн без учета их взаимодействия с торцевой стенкой. В рассмотренном выше примере после прохождения волны она равна скорости за фронтом волны ио. В дальнейшем, по мере перехода к установившемуся истечению, она уменьшается и при установившемся истечении равна нулю. Вторая составляющая (ниже будем называть ее гидродинамической) характеризует течение жидкости как несжимаемой, учитывает ее инерционность и описывает постепенный переход к установившемуся течению. В силу инерционности она непрерывна и в момент прихода волны мгновенно не изменяется. Затем она постепенно возрастает и при установившемся истечении гидродинамическая составляющая равна средней скорости установившегося истечения.
Рассмотрение акустической составляющей начнем со случая, когда волны справа и слева, распространяясь по покоящейся жидкости, встречаются в сечении дросселя. Будем условно разделять первичный фронт волны и возмущения, связанные с ростом
давления при отражении волны от торца, полагая, что последние не догоняют первичный фронт. Тогда значение иі2ак можно определить из соотношения распада разрыва по параметрам за фронтами волн:
М12ак = 7Г(М1 + и2) + 7—*—(Д_Р1 - АР2).
2 2ре
Далее, рассмотрим случай, когда волны распространяются по подвижной жидкости, причем течение близко к установившемуся (и12гидр = и12уст). При малых числах Маха можно приближенно считать, что картина взаимодействия волн с дросселем такая же, как и в случае покоящейся жидкости. Тогда данное соотношение будет справедливым, если в него подставлять избыточные (акустические) скорости и давления, равные перепадам соответствующих значений на фронтах волн. Представив эти перепады в виде разности полных и установившихся значений, получим
^12 ак = ^ (М1 ^1уст + и2 ^2уст) ^Р1уст) ( Др2 Др2уст)) •
В общем случае волны могут распространяться при неустановившемся течении. Тогда акустическая составляющая скорости определяется не только параметрами волн, которые воздействует на дроссель в данный момент (слева и справа), но и параметрами тех волн, которые воздействовали прежде, причем, чем меньше времени прошло от момента их прохождения, тем существеннее их влияние. В этом случае выделить акустические составляющие скоростей в каналах можно, заменив в последней формуле установившиеся значения на гидродинамические (см. ниже). Разность гидродинамических давлений в каналах можно выразить через потери полного напора, рассчитанные по гидродинамическим скоростям:
2 2 2 Др1гидр - Др2г„др = -р^ + р^ + Сир^^. (7)
Гидродинамические скорости в каналах пропорциональны гидродинамической скорости истечения: и^гидр = иу2гидр]г. Тогда выражение для акустической скорости примет вид
«12ак — 2 (мі _ м12гидр/і +«2 — «іггидр/г) + — Л.Р2 — Р----^^“(Сі2 + /І — /і •
Введем римановы инварианты, характеризующие акустические волны, распространяющиеся по каналам 1 и 2 в сторону дросселя: Ді+ = иі + Дрі/рс и Д2- = и — Др2/рс. Давление в каналах Дрі выразим через значения инвариантов и скоростей и*. После преобразований получим
«12ак = о ( ^1+ + ^2- — м12гидр(/і + /2)-----0 ДР (Сі2 + /| — /і ) ) • (8)
В рассмотренном выше примере непосредственно после прохождения через дроссель волны сжатия и^ак = (Д1+ + Й2-)/2 = ио, а при установлении течения и12ак ^ 0.
При определении гидродинамической составляющей скорости будем полагать, что переход от начального (акустического) режима к установившемуся определяется одной
независимой переменной, т. е. гидродинамическое течение жидкости, истекающей через дроссель, будем рассматривать как систему с одной степенью свободы. Кинетическая энергия такой системы однозначно определяется скоростью «12гидр, изменение которой обусловлено силами давления и сопротивления. С этой точки зрения гидродинамическую скорость можно было бы определить из уравнения, подобного обобщенному уравнению Бернулли (1), записанного для участка, включающего дроссель, области подтекания и струйного течения, в которое вместо полных скоростей входили бы гидродинамические. Однако следует учитывать, что при уменьшении акустической скорости по мере установления течения связанный с ней импульс не теряется, а частично переходит в импульс, связанный с гидродинамической скоростью. Этот фактор особенно существен, если площадь сечения дросселя соизмерима с площадью каналов (в рассмотренном примере — начальная акустическая скорость соизмерима с установившейся). Для описания этого фактора целесообразно ввести отдельную инерционную характеристику. Учитывая, что акустическая скорость выражается через гидродинамическую формулой (8), уравнение для последней можно записать в виде
ди12ак
^12^12гидр ^12ак"7^ ^12гидр —
ди12гидр
%^12- (9)
Здесь тою = р^12^12экв — полная эквивалентная масса жидкости, истекающей через дроссель, то12ак = pFl2 ^12ак — эквивалентная масса, характеризующая передачу импульса от акустического течения к гидродинамическому, ^12экв и Л-12ак — соответствующие эквивалентные размеры. Коэффициент ^12 рассчитывается по зависимости (4), причем при необходимости учитывается нестационарная поправка к коэффициенту трения Л [7]. Производная ди12ак/ди12гидр определяется из (8). Уравнение (9) можно рассматривать как уравнение динамики эквивалентной массы, которая разгоняется под действием силы давления (р — ^2)^12. На массу также действует сила сопротивления, которую можно представить в виде —Р^12и22гидр(/2г — /2 + С12)/2. Второе слагаемое в левой части характеризует переход импульса акустического движения в импульс гидродинамического. Скобка в правой части представляет собой разность правой и левой частей (3), т. е. в установившемся режиме уравнение (3) выполняется автоматически.
В (9) выразим гидродинамические скорости в каналах 1 и 2 через гидродинамическую скорость истечения, а давления — через римановы инварианты. С учетом введенных параметров, а также соотношений (6) и (8), уравнение (9) можно представить в виде
(/И2экв - ^ (/1 + /2 + ^Н(С12 + /| - Л2))) «12ГИдр =
с(Д1+ + Да- - «12гидр(/1 + /2)) - ^рЧС12 - /12 + /|)^ • (10)
При и12гидр < 0 коэффициент ^12 следует принять отрицательным, а в выражение (4) вместо /1 следует подставлять /2 и наоборот. В случае, если параметры, входящие в (10) (инварианты Д1+ и Д2-, коэффициент потерь £12, а также параметры Л-12ак и ^12экв) можно считать постоянными, (9) можно проинтегрировать аналитически (за
' , М1гидр , М2гидр\ А
Р1+Р—2-------\Р2 + Р—2— ) ~^12Р
нулевой момент времени примем момент прохождения волной центрального сечения дросселя). Интеграл имеет вид
^12ак ( ^12гидр ^12уст1 ^12гидр ^12уст2 ^
__ h + h ^ \^^12экв \^12гидр0 ^12уст1 ^12гидр0 ^12уст2 /
_|_ f2).^l | НС12 /i+Zf)' \^12гидр0 ^12уст1 ^12гидр ^12уст2 /
С(/1 + /2) ^ (П)
Здесь ^12гидр0 — значение гидродинамической скорости, соответствующее t = 0,
fl + h ( | 2(i?i+ + i?2-)(Cl2 - /l + /2 ) -Л Поч
у0"’2 = teTfFTTJ ( V----------------------------------------V (12)
суть корни уравнения «12ак = 0 (один из них соответствует скорости установившегося течения). Полная скорость истечения определяется из соотношения (6) с учетом (8).
Эквивалентные размеры ^12экв и hi2aK определялись из численных экспериментов (см. ниже). Расчеты показали, что их значения при постоянном направлении течения практически не зависят от его параметров. Обработка результатов 177 численных экспериментов позволяет для осесимметричного случая принять следующие аппроксимации:
^12ак = hi2 = ^12^12;
^12экв = (1 + *12)h12 + d12 (*1 (1 — /1)2 + *2 (1 — /2)“) = /10)
= d12((1 + *12 )П12 + *1 (1 — /1)2 + *2(1 — /2)“); ( )
k12 = 0,095; k1 =0, 34; k2 = 0, 39; a = 2, 3.
Здесь d,12 —диаметр дросселя. Коэффициенты *1 и *2 характеризуют приведенные массы на входе и выходе, а *12 —неравномерность течения в канале дросселя.
4. Определение средних давлений в сечении дросселя и на торцевых стенках. В некоторых задачах помимо скорости истечения необходимо также определить средние значения давления в дросселе и на торцевых стенках (последняя задача особенно актуальна, если дроссель установлен в поршне и необходимо рассчитать действующую на него силу). При установившемся истечении среднее давление в дросселе (для определенности будем рассматривать сечение в плоскости торца канала 1) можно определить из обобщенного уравнения Бернулли:
2 2 2 2 Др12Уст = АР1уст+Р^р--Р 12^СТ (^ + Cl^u) = Ар2уст+р-^-р 12^СТ (1-С12—2), (14)
где Cl—12 = Cv + C/(l-/i)£oie 11С12-2 = С12 — Ci—12 = (*V1 - /1(1 - /2) + (1 - /г)2) £oe +
+ A^12 — коэффициенты потерь на входе и выходе из дросселя. Следует отметить, что данное разделение общего коэффициента потерь на две составляющие весьма условно и требует уточнения в зависимости от параметра П12.
В случае, если по подвижной жидкости, которая истекает через дроссель в установившемся режиме, распространяется волна, давление в дросселе после ее прохождения можно представить в виде суммы установившегося давления (14), которое определяется параметрами до прохождения волны, и акустической добавки. Последняя определяется
h12
экв
t
с
1
по формуле распада разрыва, в которую вместо полных скоростей и давлений подставляются разности значений параметров после прохождения волны и установившихся параметров до ее прохождения:
2 2
АР12 = АР1уст + - Р^( 1 + Сі-і2)+
Арі Дріуст + ^Р2 ^.Р2уст и1 ^1уст (^2 ^2уст)
+ 2 + РС 2 " л и2уст и12уст л \ Ар1 Ар1уст + Ар2 Др2уст
= Др2уст + Р~Г2--------Р^—(1 ~ ^12^2) + -----------------------2--------------
и1 и1уст (и2 и2уст) /1Г,
+ рс---------2---- --------. (15)
В рамках предложенного подхода в общем случае вместо установившихся значений подставим гидродинамические. Используя (7), соотношения (14)—(15) можно обобщить на произвольный случай так:
Арі2 =
1 ( и2 и2 и2
= І ( АР1 + Р- (1 + С^Р^^ + Ар2 + р- (1 - Сі2-2)р
^Ї1\ Мігидр (^2 ^2гидр)
+ рс - .
Легко убедиться, что соотношения (14) и (15) являются частным случаем (16). Последнее можно переписать через римановы инварианты:
Д _ ^1+ + ^2- — и12гидр(/1 — І2 ) и12гидр Л /1 + /2 С12 --- 2 + С1—12
А-Р12 — рс-------------------------------------р-------- 1 —
2 г 2 V 2
При определении давлений на торцевых стенках также сначала рассмотрим частные случаи. При установившемся истечении среднее давление на торцевой стенке канала 1 может быть определено из уравнения количества движения, записанного для области канала 1, левой границей которой является сечение, в котором течение можно считать равномерным, а правой — плоскость торца дросселя. После преобразований с учетом (14) получим
Д_ _ л_ , 2/1 — /2 — 1 + ^1—12 Р^^густЛ ,лг?л
АР1СТ уСТ — АР1уСТ -\- . (1')
1 — /1 2
Если при установившемся течении жидкости со стороны канала 1 приходит волна, среднее давление на стенке можно представить в виде суммы (17) и акустической добавки, которую можно определить по формуле для отражения акустической волны от стенки, подставляя в нее (как и в (15)) разности текущих и прежних значений скорости:
Л„ _ Л„ I 2/1 — /1 — 1 + С1—12 Рм12уст/1 , ^
АР1СТ — АР1уСТ -\- + рс[и\ и 1уст)• (1®)
1 — /1 2
Приход волны со стороны канала 2 не приводит к мгновенному изменению среднего давления, поскольку расширение фронта от сечения ^12 до сечения происходит
постепенно. Обобщение (17) и (18) на произвольный случай имеет вид
Д„ _______/с ^ Ч , 2Л “ Л2 “ 1 + С1 —12 Рм12гидр/1 /1П^
Ар 1ст — рСуИ 1+ ^12гидр/1) + 1 у 2 ' (19)
И, аналогично,
ЛР2СТ = рс(ц12гиДр/2 - Д2_) + 2/2 ~ /2 ~ * ~ С12^2 рМ12г2иДР/2 . (20)
5. Сравнение результатов расчетов по феноменологической модели с численными решениями задачи в осесимметричной постановке. Для определения инерционных коэффициентов необходимо ввести критерий, по которому будет производиться сравнение с результатами численных экспериментов. Пусть в начальный момент времени жидкость была неподвижна и слева пришла волна сжатия, характеризующаяся скоростью за фронтом ио = Й1+. В качестве критерия удобно ввести безразмерное характерное время переходного процесса
_ _ с [г* и12* - и12(г) и
ТуСТ~^12./о «12* -1X12(0)
Здесь и12 (0) —начальное значение скорости истечения, равное акустической скорости ио, ию* = **и!2уст + (1 — *>12(0), к* = 0.9, и12уст — скорость установившегося течения, Ь* —корень уравнения и^^*) = и^*. При постоянном коэффициенте ^12 существует аналитическое выражение для величины туст, которое из-за громоздкости здесь не приводится.
Для тестирования модели и определения коэффициентов была разработана программа численного решения нестационарных уравнений слабосжимаемой жидкости в осесимметричной постановке. Задача решалась методом Годунова с ограничителем Чакравати—Оши. Трение на стенках не учитывалось и искусственные местные сопротивления не вводились. Для описания турбулентности во внешнем течении использовалась модель к — е (рассматривался случай больших чисел Рейнольдса). Размеры ячеек выбирались из расчета не менее 50 ячеек на радиус дросселя. Внешние границы в каналах 1 и 2 удалялись от стенок дросселя на расстояние не менее двух диаметров каналов. На границах задавались значения римановых инвариантов. Проверка адекватности результатов самих численных экспериментов заключалась в согласовании установившихся значений средних скоростей истечения со значениями, полученными по формуле (12) с учетом эмпирической зависимости для коэффициента потерь (4), а также в сопоставлении результатов, полученных при разных размерах расчетных ячеек. Погрешность численного алгоритма по отношению к значениям (12) в большинстве случаев не превышала 1-2% и лишь при П12 > 1 возрастала до 10-18%, что связано с влиянием схемной вязкости перекоса в области спрямления течения в дросселе. Для того чтобы привести модель в соответствие с условиями численных экспериментов, при определении инерционных коэффициентов по методу наименьших квадратов в выражении (4) не учитывались составляющие, характеризующие молекулярную вязкость. При существенном рассогласовании параметры, входящие в выражение для туст, определялись по откорректированным значениям коэффициента потерь, соответствующим (по формуле (12)) скорости установившегося истечения, полученной в численном расчете.
На рис. 3 показаны зависимости величины туст от геометрических параметров круглого дросселя для некоторых их сочетаний. Сплошными линиям показаны результаты
расчета по описанной модели без коррекции коэффициента потерь, пунктирными — с коррекцией (с интерполяцией между точками), а точками — результаты численных экспериментов. На рис. 4 показаны графики временных зависимостей средней скорости истечения, давления на торцевых стенках и в сечении дросселя, полученные по феноменологической модели и численным путем для одного из вариантов расчета.
Сопоставление результатов расчетов по данной модели с результатами численных экспериментов показывает, что значительные погрешности имеют место при /712 — тах(/1,/2). В этом случае течения на входе в дроссель и выходе из него формируются раздельно, т. е. рассмотрение течения через эти сопротивления как единого не вполне корректно — необходимо отдельно рассматривать каждую из этих зон и канал 12. На рис. 3, а видно, что в случае /1 — 1, /2 — 1 и П12 > 0 расчетное время установления течения начинает резко возрастать. Это связано с некорректностью предположения о неравномерности скоростей в дросселе: при шш(/1, /2) > 0, 94 вместо *12 = 0,095 следовало бы принять *12 — 0. Однако фактическая погрешность расчета скорости в этом случае ничтожна, поскольку основное увеличение скорости происходит мгновенно при прохождении волны, а величина, на которую скорость увеличивается на «гидродинамической» стадии, стремится к нулю. С другой стороны, по указанной выше причине данная модель при таких параметрах вообще является некорректной.
Рис. 3. Сопоставление характерного времени переходного процесса, рассчитанного по феноменологической модели без коррекции коэффициента потерь (сплошные линии), с коррекцией коэффициента потерь (интерполяция — пунктирные линии) и по результатам численных экспериментов (точки) для различных сочетаний параметров: а) д,\ = ^2, ио/с = 0, 00337; б) ^12/^1 = 1/6, ио/с = 0, 00337.
Также существенные отличия при использовании неоткорректированных значений коэффициента потерь наблюдаются при П12 > 1 (см. рис. 3, а), однако в этих же случаях имеет место расхождение по скоростям установившегося течения между численными расчетами и эмпирическими соотношениями (4). Вместе с тем, при использовании откорректированных значений достигнуто хорошее совпадение.
Поля скоростей при установившемся течении, полученные из численных расчетов, могут быть использованы для получения теоретических оценок используемых в модели составляющих инерционных коэффициентов на основе интегрирования кинетической энергии в соответствующих зонах. Расчеты показывают, что такие оценки дают весьма хорошее совпадение с соответствующими слагаемыми формулы (13) для составляющих, характеризующих инерцию жидкости в зоне подтекания и канале 12. Инерционный коэффициент для зоны истечения близок к полученной аппроксимации лишь
Рис.4• Результаты расчета истечения жидкости (с = 1312 м/с, р = 862 кг/м3) через тонкое отверстие (^х = d2 = 0, 3 м, dl2 = 0, 12 м, Н\2 = 0); В.2— = 0. Со стороны канала 1 последовательно приходят волны сжатия (рсК\+ = 50 • 105 Па), разрежения (рсК\+ = —50 • 105 Па) и сжатия (рсЯх+ = 0): а) средняя скорость истечения, б) средние давления в отверстии (Дрх2), на левой (Дрхст) и правой (Др2ст) стенках. Толстыми линиями показаны зависимости, полученные по описанной модели, тонкими — из численного решения уравнений гидродинамики.
при \f~fi > 0, 5. При меньших \J~fi значения, полученные путем интегрирования кинетической энергии, в 2-15 раз превышают аппроксимацию. Это объясняется тем, что в области истечения имеются зоны, скорость потока в которых, хотя и существенно, отличается от п2, но не связана напрямую с текущей скоростью истечения. К ним относятся циркуляционная зона (см. рис. 2, г) и область струйного течения на удалении более (1 — 2)с112 от дросселя. Изменение скорости в этих областях при изменении скорости истечения происходит с большой задержкой, что не позволяет механически «привязывать» кинетическую энергию жидкости в этих зонах к кинетической энергии истекающей жидкости. В численных расчетах это проявляется в медленном и весьма незначительном (1-8%) изменении скорости истечения при £ = 5 + 20£*.
Разработанная модель является основой для создания модели течения жидкости через клапаны.
Литература
1. Гиргидов А. Д. Механика жидкости и газа (гидравлика): учебник для вузов. СПб.: Изд-во СПбГПУ, 2002. 545 с.
2. Идельчик И. Е. Справочник по гидравлическим сопротивлениям. М.: Госэнерго-издат, 1960. 464 с.
3. Жуковский Н. Е. О гидравлическом ударе в водопроводных трубах. М.; Л.: Гос-техиздат, 1949. 103 с.
4. Чарный И. А. Неустановившееся движение реальной жидкости в трубах. М.; Л.: Гос. изд-во технико-теоретической литературы, 1951. 225 с.
5. Курков С. В. Метод конечных элементов в задачах динамики механизмов и приводов. СПб.: Политехника, 1991. 224 с.
6. Попов Д. Н. Динамика и регулирование гидро- и пневмосистем: учебник для вузов / Д. Н. Попов. Изд. 2-е, перераб. и доп. М.: Машиностроение, 1987. 464 с.
7. Попов Д. Н. Нестационарные гидромеханические процессы. М.: Машиностроение, 1982. 240 с.
Статья поступила в редакцию 1 октября 2010 г.
