Повна версія

Головна arrow Природознавство arrow СЕТОЧНО-ХАРАКТЕРИСТИЧНІ ЧИСЕЛЬНІ МЕТОДИ

  • Увеличить шрифт
  • Уменьшить шрифт


<<   ЗМІСТ   >>

СІТКОВИХ-ХАРАКТЕРИСТИЧНІ МЕТОД ДЛЯ НЕСТАЦІОНАРНИХ РІВНЯНЬ ГАЗОВОЇ ДИНАМИКИ

1. Дотримуючись [115], випишемо в звичайній тензорною записи диференціальні закони збереження маси, імпульсу і повної енергії для однокомпонентної суцільного середовища

що не залежать від конкретної моделі середовища і є незамкненою системою. Тут р - щільність, V - вектор відносної швидкості з фізичними компонентами і v 2 , Уз » ° kj - фізичні компоненти тензора напружень; Е - е + V 2 /2 - повна енергія одиниці маси, е = е (р, Т ) - питома (на одиницю маси) внутрішня енергія; Q - об'ємна щільність джерел (стоків) тепла; W = | w l9 w 2 , w 3 j = K r grad T -потік тепла за рахунок теплопровідності для ізотропного середовища, - коефіцієнт теплопровідності, Т - температура; в правих частинах рівнянь імпульсу і енергії включають можливі масові сили, члени, пов'язані з криволінійністю прийнятої ортогональної системи координат х х , х 2 , х 3 , а в разі, якщо система координат x lf х 2 , х 3 неінерціальна, також члени, пов'язані з переносною швидкістю руху цієї системи щодо деякої інерціальної системи координат. Зокрема, в інерціальній декартовій системі координат х, у, z при відсутності масових сил / у = f E = 0.

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

рівняння стану, що пов'язує, наприклад, тиск р, щільність р і внутрішню енергію е

умова відсутності переносу тепла за рахунок теплопровідності

У разі необхідності використовуються також калоріческой рівняння стану е = е (р, Т) і інші термодинамічні співвідношення. У цьому розділі в якості однієї з шуканих функцій вибирається питома ентальпія І = е + р / рів відповідно до цього, рівняння стану (5) використовується у вигляді

зокрема, для досконалого газу

де к = c p / c v = const - показник адіабати.

170

Вибираючи в якості шуканих величин вектор і = р 9 v v 2 , v 3 , h, запишемо тривимірні нестаціонарні рівняння течії газу (1) - (7) у вигляді (1.2.2) з гл. 1

зокрема, в разі досконалого газу (8)

Вектор правих частин f = 1 / о> / 1 * / 2 »/ з> АI» як зазначалося, залежить від вибору системи координат і при відсутності масових сил для інерційних ортогональних криволінійних систем координат з коефіцієнтами Ляме Н 1. І ,. Н * має вигляд

2. При надзвуковому просторовому обтіканні досить гладких затуплених тел кінцевих розмірів або обмежених в поперечному напрямку потоком газу, параметри якого р ж , , р ж , і т.д. є відомими функціями точки г і часу /, спостерігається картина перебігу з відійшла ударною хвилею, яка відділяє область, в якій проявляється вплив тіла ( '' обурена область "або ударний шар) від всього іншого потоку ( '' невозмущенная область") (рис. 5.1 ). Ця примикає до тіла область збуреної течії може змінювати з часом свою форму внаслідок руху її кордонів - поверхні тіла ОА ... Л7, зміни в часі параметрів набігаючого потоку або параметрів вдуваемого через поверхню тіла на деякій ділянці А г А 4 газу і т.д . Для тел порівняно простої форми в обуреної частини течії поблизу затупления є єдина область дозвукового течії, за якою розташована область надзвукової течії. Якщо в деякий момент часу t = 0 відомо положення поверхні тіла, ударної хвилі і розподіл газодинамічних функцій, то плин в будь-який момент часу /> Про повністю визначається завданням руху поверхні тіла. Для визначення цієї течії необхідно знайти рішення рівнянь (7), (9) - (11), яке задовольняє початковим умовам при t = 0 і граничним умовам на поверхні тіла і ударної хвилі, положення якої заздалегідь невідомо і підлягає визначенню разом з течією.

Для більш складних течій, наприклад, при надзвуковому обтіканні затупленого тіла, з частини якого А ь ОА 4 здійснюється інтенсивний вдув газу (в загальному випадку іншого складу, ніж в набігає потоці) в ударний шар, картина обтікання суттєво ускладнюється. Між ударною хвилею 0 3 Ci ... З 5 і тілом з'являється поверхню контактного розриву 0 2 В - В 7, яка відокремлює набігає газ від вдуваемого і, так само як ударна хвиля, що підлягає визначенню. Основна область дозвукового течії А 2 В 2 В 3 C 2 C 3 B 4 B S A S приймає складну форму, за ділянкою А 3 А 4 поблизу тіла можуть виникнути області зворотних течій, при дуже сильному вдувом можлива поява складної системи внутрішніх стрибків ущільнення, поверхню контактного розриву може замикатися на поверхні тіла і т.д.

Щоб охопити по можливості всі основні типи граничних умов, що виникають в задачах надзвукового обтікання тіл, будемо припускати, що область збуреної течії (ударний шар) складається з області I, обмеженою контактним розривом і ударною хвилею, області II між поверхнею тіла і контактним розривом, в обох цих областях існує ділянку чисто надзвукової течії і що інтенсивність інших можливих розривів істотно менше інтенсивності ударної хвилі і контактного розриву.

Для визначеності в подальшому також прийнято, що пов'язана з тілом криволинейная ортогональна система координат (з початком координат - точкою 0 всередині тіла) обрана таким чином, щоб в області інтегрування координата х 2 зростала в напрямку потоку в ударному шарі, координата х 3 відповідала різним меридіональним площинах, а координата х 2 мала напрямок, в цілому відповідне зовнішньої нормалі до поверхні тіла (рис. 5.1).

Найбільш часто буде використовуватися сферична система координат х, = s, х 2 = r, х 3 = <р, пов'язана з декартовой х, у, z співвідношеннями

Мал. 5.1

В цьому випадку Яi - г, Н 2 = 1, Н 3 = г sin s.

а) Для осесиметричних тіл і просторових конфігурацій (в тому випадку, якщо є площину симетрії течії) досить обмежитися розглядом однієї з симетричних частин течії, для якої передбачається зміна х 3 в межах між 0 і л. Тоді в площині х 3 = 0 і х 3 = п граничні умови мають вигляд

  • 147
  • 173

У розрахунках обтікання неосесиметричних тел при наявності кута ковзання р (рис. 5.1) координата 0 <д: 3 <х ^ = 2яі використовується періодичність рішення по цій координаті

Зручно і в разі (13) використовувати той факт, що рішення може бути певним чином продовжено за площині х 3 = 0 і х 3 = я, а саме

  • б) На ударній хвилі 0 3 З Х ... З 5 з деяким рівнянням х 2 х 3 )
  • (Яке, як зазначалося, підлягає визначенню) слід дотримуватися законів збереження маси, імпульсу і повної енергії, які запишемо в вигляді

- швидкість набігаючого потоку газу щодо поверхні ударної хвилі в напрямку її нормалі;

- нахили ударної хвилі до осей х % 9 х 3 відповідно; орт обраної системи координат.

Відсутні співвідношення для визначення входять в (15) невідомих р, Л, р, v 2 , в, q lc , q 3c , R c можна отримати, залучаючи рівняння стану (7), три очевидних геометричних співвідношення [44] (див. Також розд. 2 гл. І)

що дозволяють визначити нове положення ударної хвилі, а також відповідну лінійну комбінацію вихідних рівнянь (9). Для найпростішого рівняння стану (8) величина

явно виражається через в .

в) На ділянці поверхні тіла А ... А 3 Л 4 ... А 7 з відомим рівнянням * 2 = RwV, x i> * з) в відсутність вдуву граничною умовою є умова непротеканія

  • - одинична нормаль до поверхні тіла.
  • г) На ділянці А 2а поверхні тіла зі вдувом, в залежності від фізичної природи появи вдуву (витікання газів з сопла в ударний шар, руйнування теплозахисного покриття під дією великих теплових потоків і т.д.) можливі різні постановки граничних умов. При цьому слід враховувати обмеження на допустимий вид граничних умов, які випливають з розд. 5, гл. I. Зокрема, якщо механізм вдуву обумовлений досить великою величиною теплового потоку на одиницю поверхні тіла q (t, Х, ДГ 3 ), вдув здійснюється по нормалі n w до поверхні тіла, тобто V = Vn w = Vit ^ + іее ^, і

є дозвуковим (V / а <1, де а (р, h) - швидкість звуку), все термодинамічні процеси рівнозначні і відомо відповідне рівняння Клайперона-Клаузіуса для залежності тиску насичують парів від температури або ентальпії

при малих в порівнянні з ефективною теплотою руйнування X втрати на випромінювання (нагрівання поверхні тіла до температур руйнування, внутрішня і кінетична енергія пари, що утворилася входять в величину X) і при швидкості теплової хвилі (за рахунок теплопровідності), істотно меншою в порівнянні зі швидкістю хвилі руйнування

(р, - щільність матеріалу тіла в твердій фазі), то така проста модель для визначення шуканих параметрів р, р, h, v lt v 2 , v 3 (в разі, якщо D <V і R w (t, Х , х 3 ) - відома функція) включає рівняння (7), (22) і наступні співвідношення:

Останній вираз виписується на основі аналізу балансу енергії на поверхні тіла (хвилі руйнування) і саме його різні форми відрізняють ті чи інші моделі руйнування. У більш складних моделях часто враховується також нерівноважний характер процесу руйнування (в тонкому пристеночном Кнудсеновском шарі) і більш детально - процеси в твердій фазі. Для замикання системи (7), (22), (24), (25) слід залучити певну лінійну комбінацію рівнянь газової динаміки (9). У тих ситуаціях, коли порушується умова дозвукового характеру вдуву, як замикає умови следут покласти

так як '' природний "вдув не може бути надзвуковим.

Для визначення формозміни поверхні тіла, величини виносу маси і т.д. необхідно залучати закон збереження маси (23) для N w = D:

з якого можна визначати функцію R w (t, x it х 3 ) і її просторові похідні.

Замість (22), (25) можуть бути задані величини

  • (В дозвуковом випадку '' примусового "вдуву) або навіть повністю вектор і в разі надзвукового вдуву
  • д) На контактному розриві В якщо шукане рівняння його

поверхні є х 2 = R k (t, xi, x з), нормаль n k = (-qik e x l + е х 2 -

- <7з * е х,) // 1 + Як + язк> нахили до осей х х 3 , відповідно Чк = 2 / Нi) (dR k (t, x if х 3 ) / дх 1 ), q 3k = (Я 2 / Я 3 ) (ЕД до / Ех 3 ), газодинамічні параметри по обидві сторони від цього розриву Uj і иц, то граничними умовами будуть умова непротеканія

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

е) При чисельному рішенні задачі про обтікання передньої частини затупленого тіла область інтегрування обмежують зазвичай деякої замикає поверхнею { У iCiC s C ^ B 6 А 6 А 2 ). При певному її виборі (просторовий тип цієї поверхні [17,18]) на ній не потрібно постановка граничних умов. Умовою просторовості цієї поверхні є, наприклад, виконання в кожній її точці в кожен момент часу нерівності

або, в координатної записи,

і стаціонарних параметрах вдуваемого газу, також можна використовувати описану вище постановку задачі. У цьому випадку в якості початкових даних вибирають деякі досить довільні функції R? t

Rc у u °, взагалі кажучи не відповідають (9), і відшукують рішення (9), (7) при стаціонарних умовах в набігає потоці газу (37) і стаціонарних умовах вдуву. Тоді в межі при / - * 00 можливий вихід на шукане стаціонарне рішення, тобто відшукується рішення виду

Цей прийом рішення деяких стаціонарних задач називають методом встановлення і широко використовують в обчислювальній практиці.

3. Вибір системи координат в порівняно простих випадках не має принципового значення і визначається в основному міркуваннями зручності. Однак в ряді випадків область течії може мати складну геометрію з багатьма поверхнями розриву, зонами великих градієнтів і т.п. У цьому випадку або використовують складні координатні системи, що адаптуються до картини перебігу і пов'язані з поверхнями разрива- (наприклад, [116, 117] та ін.) Або з шуканим рішенням у всій області (наприклад, [118, 119] та ін.) , або будують однорідні консервативні схеми [68], що дозволяють з необхідною точністю наскрізним чином розраховувати течії при наявності особливостей всередині області інтегрування (див., наприклад, [19]). Мабуть, найбільш ефективний шлях в подібних ситуаціях - поєднання обох цих підходів з явним виділенням основних особливостей перебігу (головний ударної хвилі і т.п.) і наскрізним розрахунком розривів порівняно невеликої інтенсивності (внутрішні перегони ущільнення, хвилі розрідження і т.д. ). Не зупиняючись на детальному обговоренні цих питань, що становлять самостійний інтерес, скористаємося тут найпростішої заміною незалежних змінних типу (1.2.10), що переводить криволинейную область інтегрування в більш зручний для чисельних розрахунків прямокутний паралелепіпед з фіксованими межами (при вибоце як замикає повеохності ( 'ЗЗ 'П:

Остаточно сформулюємо розглянуту задачу. В області

знайти рішення (40), (7) і функції Л з (/, Я * (7>?>?)> задовольняють граничним умовам (13) або (14), (15) - (20), (21), ( 29) - (30), а також (22) - (25) або (24), (28) і деяким початковим умовам (35), (36).

4. Вводячи в області (42) рівномірну за координатами?, *?,? з (39) разностную сітку (рис. 5.2)

і позначаючи значення шуканих функцій в її вузлах через ь пт1к> R " mk > R пКТК »••• Для розрахунку внутрішніх вузлових точок, можна скористатися будь-який з розглянутих раніше явних різницевих схем, зокрема, недівергентной схемою I з гл. III, тобто (4.4.4). Для цього достатньо знайти відповідні власні числа і власні вектори матриць А | з.

Зокрема, для матриць А х - А е з (9), (41) маємо

Безпосередньою перевіркою можна переконатися, що А / = il ~ 'AjQ./(i = = 1, 2, 3), де Л / =; Х / (- діагональні матриці з власних значень матриць А /.

Для побудови розрахункових формул на поверхні ударної хвилі

з яких формуються замикають умови для розрахунку граничних точок. Верхній знак береться при X 2 > 0, нижній - при Х ; 2 <0.

При г? = 1 (вузлові точки т = 0, 1, ..., М, / = L , до = 0, ..., К ), з огляду на, що швидкість поширення ударної хвилі по '' рідким '* часткам в напрямку своєї нормалі більше місцевої швидкості звуку за ударною хвилею, тобто що -а < 0 <0, власні значення матриці Л 2 мають такі знаки:

і як відсутнього співвідношення до разностной апроксимації граничних умов (15) - (19) використовуємо (46) при 1 = 1. Разностная апроксимація (19) служить для визначення геометричних параметрів поверхні ударної хвилі Уестк * а підстановка (15) в

(7), (46) дає в загальному випадку нелінійну щодо систему з двох рівнянь в кожному сітковому вузлі на поверхні ударної хвилі

яка вирішується звичайним чином, тобто методом ітерацій типу методу Ньютона.

На ділянці поверхні тіла без вдуву (rj = - 1, / = 0) з урахуванням умови непротеканія (21), яке в нових змінних набуває вигляду і = 0, маємо

тобто для розрахунку таких точок необхідно скористатися співвідношеннями (46) при I = 2, 3, 4, 5 і разностной аппроксимацией умови непротека- ня (21). В результаті отримаємо в кожному сітковому вузлі, що належить цій ділянці поверхні тіла, лінійну щодо шуканого вектора u mo * систему з п'яти рівнянь, з якої нескладно отримати явні розрахункові формули для компонент цього вектора.

На ділянці поверхні тіла зі вдувом v = -H 2 R ^ t + V / l + q w + q w > > 0, отже, при дозвуковом вдувом

тобто для розрахунку цих точок необхідно скористатися одним з співвідношень (46) при / = 5 і разностной аппроксимацией граничних умов (22) - (25) або (24), (28), а при звуковому або надзвуковому вдувом - повністю задати вектор

У кожній точці контактного розриву (т? = 0, / = L / 2, т = 0, 1,... , М, до = 0, 1,..., К) необхідно визначити нове значення R ^ mk і Ю компонент векторів (і ^ д, / 2, *) ь ( u m + , I / 2, fc) n * ® кожної такої точки маємо v = v u = 0, тобто

і для їх розрахунку необхідно використовувати 8 співвідношень (46) при / = = 2, 3, 4, 5 в області I і при i = 1, 2, 3, 4 в області II, а також разностную апроксимацію граничних умов (30), (31). Разностная апроксимація (31) служить також для визначення. Наступним чисельним диференціюванням R ^ k ( або > аналогічно, (19)) можна визначити просторові Похідні Я * КТК> Ч'ктк-

Точки, що належать замикає поверхні Л Х В Х СС ^ З ^ В ь А ь А п

  • (т = М, I = 1,2 ..... L - 1, до = 0, 1,..., К) , при виконанні умови
  • (34) розраховуються, як і внутрішні вузлові точки, оскільки X! = = (Uj + а) / Н 1 > 0, А / = v x / H x > 0 (i = 2, 3, 4), XJ = (у, - а) / Н х > 0, і в співвідношення (4.4.4) увійдуть значення і на шарі t = t n = пт тільки в вузлових точках 1 М - 1, /, до I , М у I ± 1, до , i Л /, /, до ± 1 J, i М, /, А :), тобто в граничних і внутрішніх точках області інтегрування.

Розрахунок вузлових точок до = 0 або до = К при m = 1,2, ..., М - 1, / = 1, .... . . , L - 1, що належать поверхням х 3 = 0 і х 3 = xj, також проводиться у разі формулам для внутрішніх вузлових точок з урахуванням граничних умов (13) або (14). При цьому вводяться два додаткових шару вузлових точок к = - 1 н к - К + 9 параметри в яких визначаються за даними в симетричних точках к = ик = К-ь відповідно до граничних умов. Зокрема, для умови (13) маємо

У вузлових точках, що належать останньої грані паралелепіпеда (42) т = 0, / = 1,. . . , L - 1, до = 0, 1,. . . , К, при виборі системи координат (12) є особливість, пов'язана з цією системою координат. Для її усунення при розрахунку цих точок робиться формальний перехід до сферичної системі координат, поверненою щодо основної на я / 2 навколо осі у (рис. 5.1), в якій проводиться розрахунок вузлових точок m = 0, / = 1,. . . , Z, - 1, А: = 0. Параметри в точках т = 0, / = 1, ..., /. - 1, до = 1,. . . , До визначаються з урахуванням того факту, що, по суті, вони належать одній і тій же фізичній променю 0 3 (jcj = 0).

Слід звернути увагу, що в кількох найближчих у напрямку Xi до цього променю вузлових точках m = 1,2, ... 3 ) т1 до = х 2 / sinx lm ^ ** mhiX 2 i і для забезпечення апроксимації вихідних рівнянь різницевими виразами (4.4.4) необхідно виконання умови І 3 h, що не потрібно по реальній поведінці рішення в залежності від координати х 3 . Для усунення цього обмеження у напрямку х 3 в цих точках використовувалася апроксимація з другим порядком точності, що забезпечувало в цілому перший порядок точності, включаючи найближчі до променю * i = 0, лг 3 = 0 вузлові точки.

 
<<   ЗМІСТ   >>