Найпершою розрахунковою моделлю можна вважати модель основоположника пульсуючих двигунів Arhus доктора Шульц-Грунова – це знаменитий метод характеристик. Насправді дуже громіздкий і мало придатний для масового розуміння та застосування. Бажаючі можуть навіть ознайомитися з цим методом на нашому сайті, завантаживши відповідний звіт NACA, написаний цим видатним ученим. Хоча цей метод має своїх прихильників досі, він створений у 40-х роках минулого століття для ручного розрахунку на логарифмічній лінійці або, у крайньому випадку, на калькуляторі, тому мало придатний для створення програми. Принаймні про доступні для користувачів програми, засновані на методі характеристик, поки що нічого невідомо.

Те саме можна сказати і про інші методи. Деяка кількість відомих наукових праць присвячена термодинамічних моделей в рамках так званої «поршневої» аналогії руху газу в резонансній випускній трубі, коли газ розглядається як «рідкий поршень». Одним із основоположників даного методу слід визнати радянського вченого професора Є.С.Щетінкова, розробки якого належать ще до 40-х років минулого століття. До речі, аналогічний принцип був прийнятий пізніше при дослідженні пульсуючих двигунів та іншими вченими. Однак ніяких практично корисних для незалежних користувачів моделей, заснованих на аналогії «рідкого поршня», а тим більше готових програм, теж досі не з'явилося.

Відомо також про методи, засновані на 1-мірному уявленні течії газу в резонансній трубі. А ще є методи, побудовані на акустичній аналогії коливань, які виникають у резонаторі Гельмгольця. Але з усіма цими методами знову та сама проблема – звичайний користувач не може скористатися цими досягненнями і такими просунутими моделями. Тобто, користувачеві черговий раз пропонується заново повторити всі викладки того чи іншого вченого, самому скласти програму, налагодити її, повторивши всі зусилля автора моделі, і лише після цього виконати моделювання. Хто ж буде при здоровому глузді та тверезій пам'яті все це повторювати, коли і за чиї гроші? Відповідей від вчених досі не отримано...

Але найважча ситуація виникла з появою та широким впровадженням у наукову практику 3-D моделювання. Оскільки після придбання відповідної програми нічого нового чи свого винаходити вже не потрібно, деякі найбільш просунуті вчені кинулися наввипередки показувати у своїх статтях красиві кольорові картинки розподілу параметрів газу по трубі пульсуючого двигуна. Інші поспішили зосередити свої зусилля на особливостях моделювання тривимірної течії з хімічними реакціями пульсуючого горіння, якщо хороша і дорога програма це дозволяє. Треті, мабуть, у глибині наукової душі не довіряючи стандартним програмам, терміново почали самі розробляти якісь свої власні 3-D моделі пульсуючого двигуна, четверті навпаки, після не надто успішних спроб створити щось своє повернулися до стандартних 3-D програм. Правда, знову з не цілком виразною метою – чи для отримання характеристик, які вони так і не отримали, чи для того, щоб показати іншим вченим своє особливе знання просторового моделювання по чужих моделях. П'яті...

Вражаюче, але протягом 80 років тисячі вчених, маленьких і великих, досліджували пульсуючий повітряно-реактивний двигун вздовж і впоперек, написали тисячі статей, придумали сотні моделей його розрахунку і навіть видавали великих книг, але жоден з них чомусь не спромігся написати формулу для сили тяги цього двигуна. А тим часом, це одне з головних питань теорії повітряно-реактивних двигунів, саме з нього починається вивчення двигунів, а величина тяги не тільки є головною характеристикою, а й має бути метою будь-якого дослідження двигуна. Зрозуміло, якщо двигун призначений для літального апарата. Однак, зважаючи на все, це питання вченим було взагалі не цікаве, їхні роботи виконувались з якоюсь іншою метою.

Тому цілком закономірно, що всі ці глибокі дослідження об'єднує той самий результат – жодних програм моделювання пульсуючого двигуна та його характеристик, придатних для широкого практичного застосування, після таких титанічних зусиль знову не з'явилося. Користувачеві вкотре було запропоновано фактично самому купити дорогу 3-D програму, самому в ній розібратися, пройшовши відповідну підготовку, потім самому створити в цій програмі модель, налагодити її і тільки після цього, щоб напевно, моделювати свій двигун. При цьому потрібно було повторювати якісь кроки, вже виконані тим чи іншим ученим і описані у статті. У той час як повторити все це просто неможливо. Тоді виникає резонне питання – навіщо користувачеві витрачати гроші на якусь програму, час на її вивчення, створення моделі, її налагодження та моделювання двигуна, якщо у нього в гаражі є болгарка для різання листової нержавіючої сталі та зварювальний апарат?

Підсумовуючи, можна сказати, що ось уже 30 років ми спостерігаємо абсолютно парадоксальну ситуацію – вчені не можуть і, мабуть, не хочуть створити нормальну робочу модель і програму для моделювання пульсуючого повітряно-реактивного двигуна з тих чи інших причин, а користувачі не можуть і не хочуть використовувати те, що пропонують їм вчені. Більше того, жодна зі створених моделей не є працездатною через відсутність готових програм, надмірну складність та/або вартість, неможливість придбання або з якоїсь іншої причини, що перешкоджає її широкому використанню рядовим користувачем. І тоді є всі підстави стверджувати, що жодних моделей та програм для розрахунку робочого циклу пульсуючого повітряно-реактивного двигуна на сьогоднішній день фактично не створено.

І саме це визначило мету нашої роботи – створення математичної моделі та реальної програми для моделювання робочого циклу пульсуючого повітряно-реактивного двигуна саме для пересічного користувача.

Як ми це зробили?

Перше, на що ми звернули увагу, – це програми для моделювання двигунів внутрішнього згоряння. Ці програми вже стали стандартними, а розрахункові моделі, покладені в їхню основу, реалізовані як звичайні термодинамічні моделі, в яких газ у циліндрі розглядається як середовище з рівномірно розподіленими за об’ємом миттєвими параметрами. При цьому перебіг газу в трубах, прилеглих до циліндра, розглядається як одновимірне, але з урахуванням форми (діаметра) цих каналів та місцевих гідравлічних опорів. Ну, і звичайно, використовується та чи інша модель горіння, що описує закон виділення тепла в циліндрі після випаровування та займання палива.

А чому не можна зробити те саме для пульсуючого двигуна? Не тільки можна, а й потрібно, якщо ніхто цього робити не хоче. Сказано – зроблено. Ми взяли за основу термодинамічної моделі професора О.С. Щетинкова, але суттєво її переробили, особливо з урахуванням зворотної течії у трубі, коли холодне повітря надходить у трубу і потім виштовхується з неї. Ця частина моделі на даний момент заснована на моделі течії ‘рідинного’ поршня, але ми передбачили можливість переходу на 1-мірну модель течії газу в майбутньому, після повного налагодження програми, як це зроблено у відомих моделях двигунів внутрішнього згоряння.

Але найцікавіше було виявлено пізніше. Виявилося, що течію газу і повітря у впускному трубопроводі безклапанного пульсуючого повітряно-реактивного двигуна можна також представити як рух ‘рідинного’ поршня. Принаймні при правильному математичному описі поршнева аналогія дозволяє отримати ту ж фазову затримку подачі повітря в порівнянні зі зміною тиску в камері згоряння, як і будь-який інший метод, включаючи найскладніші 3-D моделі. В результаті нам удалося створити універсальну модель як клапанного, так і безклапанного пульсуючого повітряно-реактивного двигуна, вперше об'єднавши ці моделі в одній програмі.

В результаті, наша модель отримала наступні блоки компонентів.

1.Термодинамічна модель зміни тиску та температури в камері згоряння.

Якщо покласти миттєвий тиск і температуру однаковими за виділеним об'ємом камери згоряння, можна використовувати рівняння 1-го закону термодинаміки, записаного для циліндра. Звідки випливає, що:

де p – тиск;

ентальпія газу; Q – кількість тепла, що підводиться (+) та/або відводиться (–); M – маса газу в об'ємі V; Cp – теплоємність газу при постійному тиску, R – газова стала.

Розглянемо ділянку циклу пульсуючого двигуна, коли тиск у камері згоряння починає підвищуватись внаслідок підведення до газу тепла внаслідок згоряння палива. З урахуванням рівняння стану ідеального газу

після деяких перетворень з рівняння 1-го закону термодинаміки можна отримати рівняння для температури в камері згоряння у вигляді:

Де Qw – кількість тепла, що відводиться в стінки, γ – показник адіабати.

Далі, якщо диференціювати рівняння стану газу, за умови постійного об'єму камери згоряння, отримаємо:

Після перетворень обох рівнянь неважко отримати систему з 2-х диференціальних рівнянь для температури та тиску, дозволених щодо похідної, у вигляді:

де dMa/dt,dMe/dt – витрати газу з камери та повітря в камеру відповідно; qw – питомий тепловий потік у стінки; Fw – площа поверхні стін.

Розмір швидкості підведення тепла dQ/dt визначається обраною моделлю горіння. Ми прийняли наступну модель – тепловиділення має об'ємний характер, тобто відбувається рівномірно в обсязі всієї камери згоряння, при цьому швидкість тепловиділення пропорційна швидкості надходження повітря в камеру. Тобто:

де Hu – теплотворна здатність палива; η – повнота згоряння, дозволяє врахувати вплив складу паливо-повітряної суміші на кількість тепла, що виділилося; ρ – щільність газу в камері; λ – коефіцієнт надлишку повітря; L0 – стехіометричне співвідношення.

Насправді тепловиділення немає одночасно з надходженням повітря. Навпаки, в кожний момент часу впуску в камеру згоряння надходить лише мала частина повітря, яка потім перемішується з паливом, нагрівається, випаровується паливо, і тільки після цього дана порція паливо-повітряної суміші запалюється. Саме тому наша модель встановлює фазовий зсув швидкості тепловиділення щодо миттєвої масової витрати повітря. Цей фазовий зсув є час затримки займання, яке розраховується залежно від тиску та температури згідно із законом Арреніуса. Однак стосовно нашого завдання, спроба прямого обчислення часу затримки займання дає не цілком реальні результати, викликаючи нестійкість розрахунку, через що ми перейшли до відносних параметрів, але при збереженні загального характеру залежності часу від тиску та температури. В результаті отримали таку формулу:

де m – поправочний коефіцієнт, що дорівнює частці циклу двигуна, що відповідає затримці займання при початковій температурі повітря в камері згоряння T0 (попередньо прийнято m = 0,3); θ – відношення часу затримки займання при температурі T до часу за температури T0 (предварительно принято θ = 0,1). Для повноти згоряння ми також апроксимували відомі дані та отримали формулу:

де η0 – максимальний коефіцієнт повноти згоряння.

В даний час, коли ще не закінчена робота з доведення моделі, теплові втрати в стінки при вирішенні рівнянь температури і тиску тимчасово не враховуються. Натомість ми поки що просто зменшуємо величину η0, приймаючи її рівною не η0 = 0,95-0,98, а η0 = 0,70-0,75. Тим не менш, температура стінок у моделі вважається окремо за відповідними рівняннями теплового балансу з урахуванням теплоємності стінок. Так зроблено у зв'язку з виявленими проблемами збіжності розв'язання системи диференціальних рівнянь при спробі включити в ці рівняння тепловий потік у стінки (проте ми працюємо над тим, щоб знайти та усунути причину такої нестійкості).

2. Течія газу трубою.

Щоб отримати розрахункові рівняння, розглянемо перебіг з дозвуковими швидкостями та перепадами тисків і припустимо, що у всіх точках труби миттєві тиск, температуру та швидкість газу однакові. Таке припущення застосовно для порівняно коротких труб і зводить модель до так званої поршневої аналогії течії газу, тобто, до моделі ‘рідинного’ або ‘газового’ поршня. Тоді рух газу трубою розглядається як рух стовпа газу деякої маси, яка має властивості інерції. Це означає, що під дією змінного за часом перепаду тиску швидкість руху газу відставатиме від тиску – приблизно так само, як і під час руху хвиль у трубі.

Розрахункова схема методу поршневої аналогії (‘рідинного’ поршня): 1 – контрольний об'єм, 2 – конус (сопло), 3 – резонансна труба.

Таким чином, очевидно, що розв'язання задачі слід шукати у вигляді рівняння швидкості стовпа газу в трубі. При цьому важливо, що швидкість, що шукається, пов'язана з термодинамічної моделлю, що описує стан газу в контрольному об'ємі, до якого приєднана труба. Основні розрахункові рівняння, що описують газодинамічні процеси у трубі в 1-мірному наближенні – це рівняння нерозривності:

та рівняння руху (збереження імпульсу):

де t, x – час і координата за довжиною; p, ρ – тиск і щільність газу; δ – частка поздовжнього градієнта тиску, що витрачається на тертя та місцевий опір.

Тертя та гідравлічний опір може бути враховано за допомогою сумарного коефіцієнта гідравлічного опору ξ, тоді рівняння руху може бути представлене, як:

де ξ=ξ_T+ξ_c, ξ_T=λ_T L⁄D_a – коефіцієнт втрат на тертя у трубі; ξc – коефіцієнт місцевого опору у місці приєднання до об’єму. Коефіцієнти ξT і ξc при розрахунку течії в трубі можуть бути знайдені приблизно за формулами, що застосовуються для стаціонарних потоків.

Допустимо, відбувається витікання газу, що знаходиться в контрольному об'ємі під підвищеним тиском, через трубу в навколишнє середовище. Якщо рівняння руху почленно помножити на dx, то внаслідок прийнятих припущень отримаємо:

Тепер потрібно інтегрувати рівняння від центру об'єму, де приблизно можна прийняти, що газ нерухомий, до зрізу труби:

де xK – відстань від зрізу сопла до центру камери згоряння.

Після деяких перетворень можна дійти такого виду рівняння:

де

а безрозмірний параметр

є відносним об'ємом труби.

Оскільки повітря, що знаходився в трубі в початковий момент часу, швидко виштовхується газом, що має більш високу температуру, робимо припущення про те, що змішання між газом і повітрям відсутня. Тоді межа ‘повітря-газ’ рухатиметься трубою і в певний момент часу відійде від початку труби на відстань x. У трубі знаходяться гази з різними густинами, тому знайдемо деяку середню (наведену) густину газу у вигляді:

де ρx – щільність газів у трубі; ρ_0^' – щільність повітря у трубі. Якщо припустити, що за деякий проміжок часу ρ_pr пропорційна ρ, тобто: ρ_pr=rρ , де r=const, то:

Щільності газів, що входять до цього рівняння, запишемо з використанням поранення стану ідеального газу як:

де параметри з індексом ‘0’ відповідають навколишньому середовищу. За відсутності теплових втрат можна також припустити, що щільність газу об'ємом адіабатично залежить від тиску. Тоді відповідний член у рівнянні руху можна подати у вигляді:

Звідки після відповідних перетворень отримаємо остаточно при 〖p≥p〗_a і v_a≥0

Це рівняння справедливе для позитивного напрямку швидкості газу в трубі при тиску в камері більше статичного тиску на зрізі труби. Аналогічно можна отримати рівняння інших умов зміни в циклі двигуна, включаючи негативний перепад тиску і негативний (у бік камери згоряння) напрямок руху газу. Більш того, абсолютно аналогічне рівняння можна отримати і для швидкості ve газу/повітря у трубі впускної безклапанного двигуна, якщо замінити довжину L на Le і прибрати конус. До речі, знаючи швидкість, можна знайти миттєву витрату повітря через трубу, яка входить до рівнянь для розрахунку температури та тиску в камері:

У той же час, привертає увагу залежність прискорення газу від довжини труби (чим труба коротше, тим газ швидше розганяється) і від кількості повітря, що знаходиться в трубі (чим більше холодного повітря, тим газ у трубі важче і гірше розганяється). При цьому межа ‘повітря-газ’ рухається трубою зі швидкістю va і може бути визначена як:

Коли координата межі x стане менше нуля, це означатиме, що в камеру почало надходити холодне повітря. Саме цей момент використовується в моделі впускної труби безклапанного двигуна як початок впуску та нового циклу.

Тобто, модель загалом правильно враховує особливості реального процесу, принаймні вона забезпечує якісну відповідність. З кількісної точки зору порівняння є в наших статтях, наведених у ‘Бібліотеці’, звідки випливає, що розбіжність між моделлю ‘рідкого’ поршня та одновимірною моделлю у досить широкому діапазоні частот процесу теж невелика. Все це дозволило використовувати цю модель при написанні програми моделювання робочого циклу пульсуючого двигуна (проте ми вже думаємо, як у майбутньому замінити цю модель на одновимірну).

3. Клапанна система.

Для початкового етапу розробки та експлуатації програми ми прийняли найпростішу модель клапана – квазістаціонарну. Це така модель, яка є найпростішою і припускає, що частота власних коливань пелюстки клапану у багато разів більша за частоту вимушених коливань, що відповідають частоті процесу двигуна. Тобто, ця модель не враховує динамічних явищ під час руху пелюстки. Незважаючи на те, що більшість науковців упевнені, що квазістаціонарна модель – це дуже погано та неправильно, ми попередньо провели власне дослідження та порівняли її з динамічною (відповідну статтю можна прочитати у нас у ‘Бібліотеці’). Результати виявилися разючими – за сумарною цикловою витратою повітря, а це головний параметр, що визначає точність всієї нашої моделі двигуна в цілому, збіг результатів не гірший за 10-15%.

У такому разі інерцією руху пелюстки і, в першу чергу, його запізненням при зміні перепаду тиску можна знехтувати і вважати завдання квазістаціонарним. Для такого завдання, очевидно, відхилення клапана повністю визначається лише двома факторами – розподіленим зусиллям на пелюсток (перепадом тиску Δp) та жорсткістю клапана. Якщо клапан є жорстко зароблену з одного боку балку, на яку діє розподілена сила q, то залежність відхилення кінця пелюстки від цієї сили виглядає так:

де E – модуль пружності матеріалу пелюстки; Ix – момент інерції перерізу вздовж осі; l – довжина пелюстки. Враховуючи, що I_x=(δ^3 b)⁄(12,) де δ і b – товщина та ширина пелюстки, а навантаження на пелюсток від перепаду тиску q=∆pl, легко розрахувати залежність відкриття клапана (відхилення пелюстки) від перепаду тиску. Всі геометричні величини, що входять до цих формул, досить просто вважаються за всіма обраними схемами і конструкціями пелюсткових клапанів.

Прийнята схема пелюсткового клапана

Вочевидь, визначення витрати повітря через клапан необхідно також знайти гідравлічний опір клапана. Для тарілчастого клапана існує наближена формула

де y – підйом клапана; d – характерний діаметр отвору, що відкривається клапаном; b – величина перекриття отвору клапаном при повному закритті.

Однак форма пелюстки впливає на перетин, що відкривається, для проходу повітря. Тому підйом по довжині пелюстки x і площа f знаходяться шляхом розв'язування рівняння для прогину пелюстки:

Зрозуміло, ця найпростіша модель використана тільки для налагодження програми, і ми вже замінили її на динамічну, яка теж докладно досліджена нами (результати дано в тій же нашій статті). Суть динамічної моделі клапана - у вирішенні диференціального рівняння руху пелюстки, що має масу та інерцію, під дією змінного перепаду тиску:

где mx, q(x,t) - розподілена маса і навантаження на пелюсток від перепаду тисків.

Ця модель більш правильно відображає реальні властивості пелюстки, включаючи запізнення відкриття, відскок від крайніх положень (ви зможете побачити це в результатах розрахунку), і особливо закриття, коли виникає зворотне закидання газів з камери згоряння. Хоча вона теж має деякі недоліки, які ми поступово усуватимемо тим чи іншим способом у наступних версіях програми. Але ми поки що не плануємо займатися 3-D моделюванням рухів пелюстки з поперечними коливаннями. Відповідно до цього, всім любителям просторових рухів пелюсток у безповітряному просторі рекомендуємо шукати інші, більш точні онлайн програми та моделі пульсуючих двигунів.

Далі залишається тільки визначити миттєву витрату повітря через клапанну систему, що входить до рівнянь для розрахунку температури та тиску в камері:

де

коефіцієнт витрати; fe – площа, що відкривається всіма пелюстками у клапанних гратах.

У безклапанному двигуні клапанів, очевидно, немає, там впуском повітря керують не пелюсткові клапани, а труба з певним діаметром та довжиною (див. вище).

4. Інтегрування рівнянь

В результаті всіх дій виходить система диференціальних рівнянь, що описує робочий цикл пульсуючого двигуна:

У безклапанному двигуні до системи додається ще одне рівняння виду:

Рішення системи виконано шляхом Рунге-Кутта 2-го порядку (модифікований спосіб Ейлера). Це 2-кроковий метод, на першому кроці (половина кроку за часом Δt/2) обчислюється проміжне значення функції (T, p, v):

Далі розрахунок ведеться спочатку з новим аргументом функції, після чого на 2-му етапі обчислюється уточнене значення функції

Ми перевірили точність інтегрування в циклі за допомогою звичайного 1-крокового методу Ейлера. З'ясувалося, що 2-кроковий метод дає лише незначні поліпшення точності інтегрування параметрів циклу, трохи більше 1-2%. Це говорить про те, що використовувати методи вищого порядку для нашого завдання просто не має сенсу.

Разом з параметрами газу інтегруються маса повітря та газу M, координати межі ‘повітря-газ’ у трубі x та імпульс, що створюється трубою (у безклапанного двигуна – двома трубами). Крім того, в моделі розраховуються й інші допоміжні параметри, у тому числі, показник адіабати та постійна газова продуктів згоряння від температури і складу суміші.

5. Розрахунок основних параметрів двигуна

Розрахунок основних параметрів здійснюється на самому початку впуску: • у клапанному двигуні – при зміні перепаду тиску в камері з позитивного (вище атмосферного) на негативне, оскільки перепад тиску визначає впуск, • у безклапанному двигуні – на початку надходження повітря в камеру (координата межі повітря-газ у впускній трубі стає менше нуля).

Логіка алгоритму визначає ці моменти з високою точністю, і перше, що потрібно, визначити час циклу tc як час початку впуску від аналогічного попереднього моменту часу. Це принципово важлива особливість моделі – цикл ніяк не задається, а визначається процесами, які в ньому відбуваються та описані відповідними рівняннями.

Далі, після завершення циклу слід обчислення параметрів циклу. Це, зокрема:

частота циклів:

витрати повітря та палива:

тяга двигуна та питома витрата палива:

а також безрозмірні потяги, витрата повітря та частота циклів (число Струхаля):

У цьому розрахунок циклу закінчується і починається розрахунок наступного робочого циклу.

Якщо різниця між наступним та попереднім циклами зменшується, то через 5-7 циклів настає або режим автоколивань (двигун працює), або згасання коливань (двигун не працює), причому різниця в тязі між циклами автоколивань обчислюється як похибка моделі - pазвичай це менше 1-2%. В результаті за режимом автоколивань можна визначити працездатність двигуна обраної геометрії, а далі, змінюючи геометричні розміри, можна підібрати оптимальні по тязі і витраті палива, тобто створити свій власний пульсуючий повітряно-реактивний двигун і перевірити його роботу до/без різання та зварювання. На сьогоднішній день переважна більшість моделей, відомих з наукових праць, особливо 3D моделі з красивими кольоровими картинками, з такою роботою не справляється.

Більш детально деякі основні положення моделі викладені у наших статтях, які представлені на сторінці ‘Бібліотека’. Ну, і в міру вдосконалення моделі ми викладатимемо тут усі зміни, включаючи посилання на наші статті чи інші джерела.