Диплом. Пример++димплома. Кваліфікаційна робота пояснювальна записка
 Скачать 1.76 Mb.
|
|
Дискретне рівняння Больцмана Щоб отримати можливість моделювати динаміку безперервного рівняння Больцмана на комп'ютері, нам необхідно його дискретизувати. Для цього спо- чатку введемо рівномірну сітку просторових координат, крок сітки нехай буде однаковим по всіх осях. Поведінку рідини ми будемо визначати саме в вузлах сітки. Фактично, ми дозволяємо молекулам перебувати тільки в певних просто- рових вузлах. Крім того, ми повинні дискретизувати час, ми будемо визначати стан рідини в рівновіддалені один від одного моменти часу. Крім того, ми до- зволимо молекулам мати тільки певні значення швидкості такі, щоб за крок за часом вони встигали перейти в який-небудь сусідній вузол. Ці дозволені на- прямки будуть однаковими для всіх просторових вузлів. Очевидно, що по діа- гональних напрямках швидкості частинок будуть більше, ніж за недіагональні. Інтуїтивно можна зробити висновок, що при нескінченно малому кроці часу і кроці просторової решітки ця дискретна система перейде в звичайне рів- няння Больцмана, яке, в свою чергу, переходить до рівняння Нав'є-Стокса в ма- кроскопічному межі. У подальшому викладі всюди передбачається, що система одиниць така, що крок решітки є одиницею довжини, а крок за часом - одиницею часу. Для простоти нижче будемо припускати, що зовнішні сили відсутні. Ми пронумеруємо дозволені напрямки швидкості від 1 доз індексом i . Якщо тепер масу частинок, що пролітають з цього вузла в напрямку i за крок часу по- значити як i f : 26 , 1 ( , ) e i i i q i i f r v t f r f f t (2.7) Тут ми врахували, що крок за часом дорівнює одиниці, і замінили всі dt з формули (1.4) на одиницю. eq i f позначає дискретну рівноважну щільність роз- поділу, що залежить від макроскопічних маси і швидкості в даному вузлі. Мине вказали, з якого саме вузла будемо використовувати eq i f - з , i r v t в момент часу 1 t або з r в момент часу Для обчислювальної схеми зручніше виявляється використовувати вузол , i r v t в момент часу 1 t . Тоді рівняння вище можна розкласти на дві складо- ві, поширення (streaming step) і зіткнення (collision step). Поширення (streaming step): , 1 ( , ) e i i i q i i f r v t f r f f t Зіткнення (collision step): , ( , ) i i eq i i f r t f r t f f (2.8) Для дискретизованих напрямків швидкості маса і макроскопічна швид- кість в кожному вузлі будуть розраховуватися як i i i i i f u f v (2.9) Тут і далі ми пишемо всюди щільність замість маси, оскільки при одини- чному просторовому кроці решітки на кожен вузол доводиться одиничний об- сяг, і маса і щільність співпадуть за значенням. 27 Ми покажемо, що рівноважні функції розподілу залежать від маси в вуз- лах і макроскопічної швидкості. Тому після streaming step необхідно перераху- вати маси і швидкості в кожному вузлі, перерахувати рівноважні функції роз- поділу, і потім здійснювати collision. Таким чином, накожному кроці обчислювальної схеми необхідно по- ширити» частинки, тобто перемістити частки, що летіли з вузла r в напрямку в вузол i r v (зробити це для всіх частинок і напрямків). Після цього необхідно перерахувати маси, швидкості, рівноважні функції розподілу. Нарешті, необхі- дно «зіштовхнути» частинки, які прилетіли в даний вузол - тобто перерозподі- лити частки за напрямками. Проілюструємо обчислювальну схему на прикладі двомірної системи. Дискретизація на просторові вузли і зв'язки між ними (тобто дозволені напрям- ки швидкості) зображені на рисунку нижче. Просторові вузли позначені кружечками, зв'язку між ними - тонкими лініями. Рисунок 2.3 – Обчислювальну схему на прикладі двомірної системи На наступному рисунку зображено одна ітерація пари «streaming / collision». Кольорові стрілки зображують потоки літаючих молекул. Інтенсив- ність кольору кодує масу молекул, що летять в даному потоці, довжина стрілок приблизно відповідає шляху, який проходить потік за крок по часу (лише приб- лизно, оскільки стрілки повинні йти від центру вузла до центру вузла). 28 Рисунок 2.4 – Ітерація пари «streaming / collision» Просторові решітки В LBM під гратами (lattice) розуміється набір дозволених векторів швид- кості (однаковий для кожного просторового вузла). Це узгоджується зі стандар- тним математичним визначенням про решітку як про сутність, за допомогою якої шляхом паралельних перенесень можна отримати просторову сітку (grid). У LBM будь-яка решітка повинна містити нульовий вектор з вузла в себе самого - він описує частинки, які нікуди не летять з цього вузла. У LBM решіт- ки зазвичай позначаються абревіатурою DnQm, де n - розмірність простору, m- число векторів в решітці. Наприклад, D2Q9, D3Q19 і т.д. У двовимірному просторі LBM решітка може складатися, наприклад, з 5 векторів (2 вертикальних, 2 горизонтальних і нульовий вектор з вузла в себе самого, а може з 9 векторів, як на рисунку вище (2 вертикальних, 2 горизонта- льних, 4 діагональних, один нульовий). Це решітки D2Q5 і D2Q9, відповідно. Очевидними факторами для вибору решіток є – точність моделювання (інтуїтивно зрозуміло, що чим більше векторів в решітці, тим точніше моделю- вання) та обчислювальні витрати (розрахунок на решітці D2Q5 буде швидше розрахунку на D2Q9). Як не дивно, цене найважливіші фактори. Найважливіші фактори - це відтворюваність рівнянь Нав'є-Стокса і симетрія деяких тензорів, побудованих на векторах решітки. Зазвичай використовуються решітки D2Q9, D3Q15, D3Q19. Решітки D2Q9 і D3Q19 зображені нижче. Базисні вектора решітки зазвичай позначають- 29 ся як i е або с (вони співпадуть зі швидкостями i v , введеними раніше при оди- ничному кроці за часом. Нижче ми будемо використовувати позначення i е Рисунок 2.5 – Решітки D2Q9 і D3Q19 Випишемо базові вектора для D2Q9: (0,0) 0 (1,0),(0,1),( 1,0),(0, 1) 1, 2,3, 4 (1,1),( 1,1),( 1, 1),(1, 1) 5,6,7,8 i і е і і (2.10) Випишемо базові вектора для D3Q19: (0,0,0) 0 ( 1,0,0),(0, 1,0),(0,0, 1) 1, 2....5,6 ( 1, 1,0),( 1,0, 1),(0, 1, 1) 7,8...17,18 i і е і і (2.11) Рівноважні функції розподілу У безперервному випадку рівноважна функція розподілу (розподіл Макс- велла-Больцмана) дорівнює: 30 2 ( ) 2 /2 (2 ) v u eq RT D f e RT (2.12) де R - універсальна газова постійна; Т - температура D - розмірність простору v - вектор швидкості. Тут молярна маса газу приймається рівною одиниці (вона для нас не важ- лива - важливо тільки макроскопічна щільність). Можна вважати це зміною одиниці маси в нашому моделюванні. Крім того, функція нормується на лока- льну макроскопічну щільність, а не на одиницю. Зауважимо також, що зазвичай від v не віднімають макроскопічну швидкість газу u . Це відбувається тому, що зазвичай розподіл досліджують в разі нерухомого газу, нам же необхідно перейти в локальну інерційну систему відліку, що рухається з поточною швидкіс- тю газу в даній точці в даний момент часу, щоб використовувати розподіл Мак- свелла-Больцмана. Віднімання u і є таким переходом. Припустимо, що в даній точці простору розподіл молекул за швидкістю виявилося рівноважним. Цей розподіл залежить від макроскопічної маси ρ і швидкості u. З іншого боку, ми можемо обчислити ρ і u по функції розподілу. Основною вимогою для дискретних рівноважних функцій розподілу є відтворення рівняння Нав'є-Стокса в межах нескінченно малих кроків за часом і кроків решітки. Це еквівалентно тому, що якщо при заданих ρ і u ми спробуємо їх знову розрахувати по рівноважним функцій розподілу в безперервному і дискретно випадку, результати співпадуть. eq eq eq eq i i i i i f dv f u f vdv f v (2.13) Для однозначного визначення рівноважної дискретної функції розподілу 31 нам необхідно буде також включити аналогічну вимогу рівності макроскопіч- них теплових енергій в даній точці, ми його опустимо для стислості. Якщо підсумувати за напрямками швидкості рівняння для collision step, то можна показати, що collision step не змінює макроскопічні маси і швидкості молекул, що знаходяться в вузлі. Виявляється, якщо просто підставити дискретні вектори швидкості в без- перервну рівноважну функцію розподілу, рівності не виконуватимуться. Виявляється також, що ми можемо змусити рівності виконуватися, якщо розкладемо розподіл Максвелла-Больцмана вряд Тейлора до другого порядку макроскопічної швидкості. Це відповідає тому, що величина u RT дуже мала. Це обмеження завжди має виконуватися в процесі моделювання у всіх вузлах. Але навіть розкладання вряд Тейлора недостатньо. Ми також повинні ввести в дискретизованої функції eq i f спеціально підібрані множники i w . Роз- рахунок відбувається через тензори, аж до четвертого рангу, побудовані навек- торах решітки, остаточно отримаємо 2 2 2 1 2 2 eq i i i e u e u u f w RT RT RT (2.14) Але тут криється підступ: ми всюди припускаємо, що крок решітки і часу є одиницями довжини і часу, відповідно. Тому мине можемо взяти значення R з СІ, а температура тут не дорівнює температурі рідини яка моделюється в Ке- львінах. Щоб визначити їх значення, зазначимо таке. припустимо, що в рідині є обурення - тобто в деяких вузлах є надлишкова маса. Ця маса почне «розтікати- ся» по простору, до того ж в крайньому правому вузлі в області обурення вона буде рухатися в тому числі і по напрямку (1, 0, 0) в 3D або (1, 0) в 2D. За оди- ницю часу молекули вздовж цих напрямків пройдуть одиницю довжини, тобто 32 їх швидкість буде дорівнює одиниці. Це означає, що швидкість звуку як швид- кість розповсюдження збурень в нашій системі теж дорівнює одиниці. У результаті 2 2 9 3 1 3 ( ) 2 2 eq i i i i f w e u e u u (2.15) Зведемо в таблицю 2.2 всі переваги та недоліки розглянутих методів Таблиця 2.2 – Порівняння методів моделювання Назва методу Переваги Недоліки Рівняння Нав'є- Стокса Дозволяє отримати параметри швидкості та роз- поділу тиску в зоні дослі- дження з високою точністю Потребує великої кількості обчислювальних ресурсів та відповідно знач- ного часу для отримання розв’язку навіть для достат- ньо простих задач. Кліткові автомати Кліткові автомати прості у використанні для моделювання складних про- цесів та об’єктів і викорис- товуються у різноманітних галузях науки - Неізотропний терм ад- векції - Порушення інваріантно- сті Галілея - Побічні інваріанти - Шум - Тиск залежить явно від швидкості Гратчасті рівняння Больцмана В даній моделі вирі- шено такі проблеми клітко- вих автоматів, як шум, зале- жність тиску явно від швид- кості, порушення інваріант- ності Галілея. Цей метод є найбільш популярним і роз- винутим. LBM може стати не- стабільним в деяких системах при високих числах Рейнольдса У результаті аналізу описаних методів моделювання в задачах гідродина- міки запропоновано для вирішення поставленої задачі використовувати метод моделювання рідин на основі гратчастої моделі Больцмана. 33 3 ПРОГРАМНА РЕАЛІЗАЦІЯ ТА РЕЗУЛЬТАТИ ОБЧИСЛЮВА- ЛЬНОГО ЕКСПЕРИМЕНТУ 3.1 Опис програми В данному розділі, показано, як створювалась модель руху човна, з ви- користанням решітчатої моделі Больцмана та бібліотеки OpenLB. Створення вікна програми на Linux, та графічне відображення всіх обєктів та процесів, ви- конане з допомогою таких бібліотек, як OpenGL, та QT. Весь процесс включає в себе такі основні кроки: 1 крок. Ініціалізація. На цьому кроці встановлюються перетворення між фізичними та решітковими блоками. Він також визначає, де зберігаються дані моделювання і який тип решітки використовується. 2 крок. Підготовка геометрії. Геометрія підключається з іншого файлу файл stl), або з визначення функцій індикатора. Потім сітка створюється і іні- ціалізується на основі заданої геометрії. 3 крок. Підготовка решітки. На цьому кроці у відповідності з номерами матеріалу геометрії встановлюється динаміка решітки. Цей крок характеризує модель зіткнення і поведінку кордону вибраної моделі. 4 крок. Головний цикл з таймером. Спочатку таймер ініціалізується. Потім починається відлік, після чого цикл над усіма тимчасовими кроками iT починає моделювання, під час якого функції setBoundaryValues, collideAndStream і getResults (описані далі кроки 5, 6 і 7) викликаються повторно, поки моделювання не буде сходитись або не буде досягнуто максимум іте- рацій. 5 крок. Визначення початкових і граничних умов. Перша з трьох важ- ливих функцій, які викликаються під час циклу, це setBoundaryValues. Ця фун- кція задає повільно зростаючу граничну умову припливу води. Оскільки грани- ця залежить від часу, це відбувається в основному циклі. Межі можуть залиша- 34 тися незмінними протягом всього моделювання, і функція неповинна робити нічого після першої ітерації. 6 крок. Зіткнення і виконання потоку. Друга функція collideAndStream викликається в кожному кроці ітерації, щоб виконати колізію і потоковий крок.. 7 крок. Обчислення та виведення результатів. В кінці кожного кроку ітерації викликається функція getResults, яка створює консольний вихід, .vti файли результатів у певні кроки виконання програми. 8 крок. Графічне відображення руху човна в річці. За допомогою гра- фічних бібліотек С+ таких як QT, та OpenGL, будується графічний інтерфейс та малюються основні складові програми. Малюється човен, течія, земля, пові- тря. З використанням результатів моделювання записаних в файли графічно ві- дображається рух човна проти течії і зміна руху води відносно човна. Далі будуть описанні детальніше всі кроки написання програми з пояс- ненням. Розглянемо початок програми. Щоб виконати моделювання та отрима- ти деякі результати, завантажимо та розпакуємо OpenLB у системі Linux. OpenLB не має інших залежностей і вимагає лише компілятора С+ по- ряд з реалізацією OpenMPI. Обидва продукти можуть бути встановлені на си- стемі Ubuntu таким чином Створимо файл river.cpp. Потім створюємо виконуваний файл, скомпілю- ємо програму за допомогою команди make. Нарешті, запусимо моделювання на ./river і спостерігаємо за висновком терміналу. Кілька рядків для програм з OpenLB будуть завжди однакові, див. Лістинг 1. 35 Лі- стинг 1. river.cpp Початок програми на OpenLB. Приведемо деякі пояснення до коду Рядок 1: Файл заголовка olb2D.h включає визначення для всього коду, присутнього у випуску. Рядок 3: Більшість коду OpenLB залежить від параметрів шаблону. Тому його не можна заздалегідь скомпілювати, і його потрібно інтегрувати як є у ваші програми через файл olb2D.hh. Рядок 11: Весь код OpenLB міститься в просторі імен olb. Дескриптори мають власний простір імен і визначають компонування решітки, наприклад. D2Q9 або D3Q19. Рядок 14: Вибір двовимірної арифметики з плаваючою точкою. Рядок 15: Вибір дескриптора решітки. Дескриптори решітки вказують не тільки, яку решітку ви збираєтеся використовувати (для 2D моделювання, по- точний випуск OpenLB не дає вам вибору, окрім D2Q9). Наступний код представляє короткий огляд структури програми OpenLB, див. Лістинг 2. Тут показані основні функції, які будуть реалізовані далі в цьо- му розділі. 36 37 Лістинг 2. Основні функції для програми на OpenLB. OpenLB має інтерфейс для даних на основі геометрії stl і генерує повністю автоматизовану регулярну сітку. Будуємо наш обчислювальний домену простих кроків: 1. Крок: Створення екземпляра Indicator2D. Читаємо файл. 38 2. Крок: Побудуємо SquareGeometry2D. Під час будівництва квадрати бу- дуть автоматично видалятися, скорочуватися і зважуватися для хорошого балансу навантаження. 3 Крок: Побудуємо LoadBalancer, який призначає квадратам потоки. 4. Крок: Побудуємо SuperGeometry2D, який пов'язує номери матеріалів з вокселями. 5. Крок: Встановимо номери матеріалам, щоб визначити динаміку для рі- дини та межі для моделювання. 6. Крок: Побудуємо SuperLattice2D для виконання алгоритму потоку та зіткнення. 39 Лістинг 3. Створення геометрії на основі STL. Всі шість кроків подають- ся коротко як вихідний код. Створення stl файлу виконуємо наступним чином. Геометрію створюємо за допомогою такого інструменту CAD, як FreeCAD. По-перше, 2D малюнок створюємо на обраній площині (площина xy), ви- користовуючи квадрати і багатокутники. Кілька таких 2D-об'єктів обєднуємо за допомогою операцій, таких як об'єднання, перетин, тощо для отримання цільо- вої геометрії. В середовищі FreeCAD створюємо такі основні деталі нашої моделі, як - Земля - Вода - Повітря - Човен Зберігаємо результат в файл geometry.stl. Розглянемо детальніше налаштування номерів матеріалу. OpenLB має загальне правило для представлення геометрії. Кожному осередку присвоюєть- ся конкретне число, яке називається числом матеріалу, визначаючи, чи лежить ця клітина на кордоні або в області рідини, та чи є вона зайвою в обчисленнях. Рисунок 6 ілюструє цена нашому прикладі потоку. Різні кроки зіткнення і потоку на кордоні і рідині визначаються по відношенню до номера матеріалу. Пе- ревагою використання чисел матеріалу в моделюванні потоку є автоматичне визначення напрямків рідини на граничних вузлах, оскільки цене завжди практично вручну, якщо номер матеріалу складної геометрії отримано з файлу stl. 40 Встановлюємо чисела матеріалів за допомогою однієї з функцій перей- менування в SuperGeometry2D. Розглянемо код // replace one material number another material number void rename(int fromM, int toM); // replace one material that fulfills an indicator functor condition with another void rename(int fromM, int toM, IndicatorF2D Лістинг 4. Функція перейменування для встановлення номерів матеріалу. Рисунок 3.6 – Звязок вузлів решітки геометрії з номерами матеріалів. Номер матеріалу 0, 1, 2 ,3 ,4 і 5 відповідають геометрії, рідини, межі повернення, припливу, відтоку і перешкоді, (1 = рідина, 2 = кордон без ковзання, 3 = межа швидкості, 4 = межа постійного тиску, 5 = вигнутий кордон, 0 = нічого). |