1.4.9. Чисельні методи розв’язання рівнянь Нав’є-Стокса нестисливої рідини
Рівняння нерозривності:
.
Рівняння перенесення імпульсу
, . (1.121)
Рівняння перенесення тепла
, . (1.122)
Рівняння Нав’є-Стокса для стисливої рідини утворюють змішану систему гіперболічно-параболічних рівнянь, а для нестисливої рідини еліптично-параболічних. Тому в цих двох випадках використовуються різні чисельні методи розв’язання.
Для розв’язання рівнянь стисливої рідини широко використовуються явні методи розрахунку: методи Дюфорта-Франкела, Лакса-Вендроффа, Мак-Кормака та ін. Використання цих методів для розрахунку нестисливих течій не можливе через обмеження на крок за часом, що накладається умовою Куранта-Фрідріхса-Леві
.
Нестисливість течії приводить до додаткових обчислювальних труднощів: до рівняннь нерозривності входять лише компоненти швидкості. Тому в даному випадку немає прямого зв'язку з тиском, який у разі стисливих течій здійснюється через густину. Для розрахунку тиску потрібно розв’язувати еліптичне рівняння Пуассона разом із рівнянням імпульсів.
Особливості дискретизації рівнянь Нав’є-Стокса
Найчастіше трапляються рівняння перенесення субстанцій, які для нестисливої рідини мають вигляд узагальненого рівняння перенесення:
. (*)
Форма запису рівнянь. Для скінченно-різницевої дискретизації можна використовувати як дивергентну форму (*), так і недивергентну форму
,
проте досвід чисельних розрахунків (наприклад, для тестової задачі течії у квадратній каверні) показав, що дивергентна форма запису (*) приводить до точніших результатів.
Замість параметризації вузлів сітки використовується одновимірна згортка тривимірних координат з використанням позначень, які відповідають сторонам світу (див. рис. 1.10, табл. 1.1).
Система скінченно-різницевих рівнянь (п'ятиточковий шаблон для двовимірного рівняння) записується у вигляді
Проінтегруємо узагальнене рівняння перенесення по елементарному об'єму, грані якого проходять через середини відрізків, що з’єднують даний вузол із сусідніми (див. рис. 1.11).
. (1.124)
Апроксимації поверхневих інтегралів: – сума по всіх гранях об'єму.
Рисунок 1.11 – Контрольний об’єм
Розглянемо грань e і інтеграл
.
1) Формула прямокутників
– 2-й порядок апроксимації (найпростіша).
2) Формула трапецій
– 2-й порядок апроксимації.
3) Формула Сімпсона
– 4-й порядок апроксимації.
Апроксимації об'ємних інтегралів
2-й порядок апроксимації.
Застосувавши апроксимації за методом прямокутників, одержимо
.
Дифузійні члени
– 2-й порядок апроксимації.
Апроксимації конвективних членів:
Інтерполяція по потоку (Upwind Differencing Scheme)
2-й порядок апроксимації, але при розв’язанні виникають осциляції через порушення умови монотонності і діагонального переважання при великих значеннях сіткового числа Рейнольдса
.
Квадратична інтерполяція по потоку (QUICK):
(1.127)
-
3-й порядок апроксимації, але умови монотонності і діагонального переважання при великих також не виконуються. Через це виникають проблеми при ітераційному розв’язанні відповідної системи лінійних рівнянь (повільна збіжність), а також нефізичні осциляції.
Крім того, одновимірний розрахунковий шаблон перестає бути 3-точковим, що не дозволяє використовувати метод прогону при розв’язанні.
Уникнути цих труднощів допомагає метод відкладеної корекції (Deferred Correction). Скінченно-різницева апроксимація високого порядку подається у вигляді суми противопотокової схеми і залишку, який виноситься у джерельний член і обчислюється з попередньої ітерації
. (1.128)
Приклад. Розглянемо рівняння дифузії
,
.
Схема UDS 1-го порядку точності:
, .
Схема QUICK 1-го порядку точності:
,
,
останній вираз можна записати у вигляді
(1.129)
Метод відкладеної корекції реалізується у вигляді ітераційного процесу
(1.130)
на кожній ітерації використовується метод прогону, рахунок за яким стійкий.
Ітераційні алгоритми сумісного розв’язання рівнянь перенесення імпульсу і нерозривності: методи проекції
(1.131)
У сім’ї методів корекції тиску проблема вирішується таким чином. Спочатку знаходиться проміжне поле швидкості з рівняння імпульсів із «старим» тиском
.
Поле швидкості не задовольняє рівняння нерозривності. Тому на наступному етапі шукаються поправки до поля швидкості і тиск для того, щоб задовольнити це рівняння: .
Підставляючи в (1.131) і в рівняння нерозривності, одержимо рівняння для поправок швидкості і тиску
.
Нехтуючи в'язкими і конвективними членами, одержимо
. (1.132)
Інтегруючи це рівняння на інтервалі , одержимо
, (1.133)
звідки, з урахуванням виразу для дивергенції поправки швидкості одержимо рівняння Пуассона для тиску, розв’язуючи яке знайдемо шукані поправки до швидкості з (1.133).
.
З (1.133) видно, що поправка до швидкості є безвихоровою, оскільки є градієнтом деякого потенціалу, який пропорційний до тиску. Методи подібного типу, в яких на першому кроці обчислюється поле швидкості, що не задовольняє рівняння нерозривності, а потім з нього виділяється соленоїдна складова («соленоїдна проекція»), називаються методами проекції. Наближенню (1.132) можна додати цікаву фізичну інтерпретацію. Рівняння (1.132) збігаються з рівняннями теорії імпульсного руху рідини (наприклад, у результаті удару). Такі рівняння випливають з повних рівнянь Нав’є-Стокса в межах дуже малої тимчасової тривалості удару, тобто коли
,
де величина називається імпульсом тиску.
Таким чином, наближення (1.132) відповідає короткочасному гідравлічному удару, який «повертає рідину в нестисливий стан».
Метод SIMPLE
Розглянемо реалізацію цієї ідеї в популярному методі SIMPLE (Semi-Implicit Method for Pressure-Linked Equation)
Скінченно-різницеве рівняння імпульсів, яке необхідно розв’язати на n+1 кроці за часом, має вигляд
. (1.134)
Рівняння розв'язується ітераціями, які називатимемо зовнішніми ітераціями (на відміну від внутрішніх ітерацій при розв’язанні системи лінійних рівнянь).
На першому кроці m-ї зовнішньої ітерації розв'язується рівняння для проміжного поля швидкості
,
потім шукається така поправка до нього, яка задовольняла б скінченно-різницеве рівняння нерозривності, тобто
. (1.135)
Рівняння для відповідних поправок швидкості і тиску випливає з (1.134):
. (1.136)
У методі SIMPLE (Patankar, Spalding, 1972) сумою в лівій частині нехтує, після чого виходить наближений розв’язок
, (1.137)
підставивши який в (1.135), одержимо рівняння Пуассона для поправки тиску
. (1.138)
Як показали чисельні експерименти, обчислену звідси поправку тиску необхідно піддати нижній релаксації так, що підсумкові вирази для тиску і швидкості мають вигляд
де .(1.139)
Рекомендоване значення становить близько 0.8 (Patankar, 1980). Рівняння для проміжної швидкості також рекомендується вирішувати за допомогою оригінального варіанта методу нижньої релаксації, запропонованого Патанкаром, в якому замість рівняння (1.134) розв’язується
,(1.140)
де – коефіцієнт нижньої релаксації, оптимальне значення якого становить (Ferziger,Peric, 1996)
. (1.141)
Фактично, рівняння (1.140) можна одержати з (1.134) додаванням в ліву частину похідної за псевдочасом, причому крок за псевдочасом дорівнює
. (1.142)
Рівняння (1.140) при цьому розв'язується неявним методом до встановлення, що підвищує стійкість ітераційного процесу.
Метод SIMPLEС
«Слабким місцем» наведеного вище рівняння для поправки тиску (1.138) є нехтування сумою в лівій частині (1.136). З погляду розглянутої вище аналогії дії поправки тиску з гідравлічним ударом таке нехтування означає, що швидкість від дії імпульсу поправки тиску змінюється тільки в точці P, тоді як для сусідніх точок ми вважаємо . В той самий час гідравлічний удар діє на весь об'єм рідини відразу, тому швидкість у сусідніх вузлах повинна змінюватися на величину такого ж порядку. Виходячи з цих міркувань, можна уточнити наближення (1.137), взявши в (1.136) . Тоді замість (1.137) одержимо
. (1.143)
Рівняння Пуассона для поправки тиску (5) набуває вигляду
. (1.144)
Даний метод відомий як SIMPLEC (SIMPLE Corrected) (Van Doormal, Raithby, 1984). Він має швидшу збіжність при розв’язанні стаціонарних задач. З порівняння (1.141),(1.142) із співвідношеннями (1.137) – (1.139) методу SIMPLE видно, що метод SIMPLEC еквівалентний методу SIMPLE з коефіцієнтом релаксації тиску, що дорівнює
. (1.145)
Неважко перевірити, що для всіх консервативних схем повинне виконуватися співвідношення на коефіцієнти різницевої схеми
, (1.146)
де – внесок від нестаціонарного члена. Для повністю неявної схеми, наприклад, , тоді рівняння для поправки тиску в методі SIMPLEC набуває такого вигляду, що збігається з рівнянням (1.136)
.
Для стаціонарних задач і необхідно розв’язувати рівняння для швидкості з нижньою релаксацією (1.140). У цьому випадку
,
і вираз (1.143) збігається з наведеним вище співвідношенням (1.141), що можна розглядати як його обгрунтування.