Размещено на http://www.allbest.ru/

НАЦІОНАЛЬНИЙ ТЕХНІЧНИЙ УНІВЕРСИСТЕТ УКРАЇНИ “КИЇВСЬКИЙ ПОЛІТЕХНІЧНИЙ ІНСТИТУТ”

Лук'яненко Святослав Олексійович

УДК 519.6

АДАПТИВНІ МЕТОДИ ЧИСЕЛЬНОГО МОДЕЛЮВАННЯ ВИСОКОГРАДІЄНТНИХ ПРОЦЕСІВ В ОБ'ЄКТАХ З РОЗПОДІЛЕНИМИ ПАРАМЕТРАМИ

Спеціальність 01.05.02 - Математичне моделювання та обчислювальні методи

Автореферат

дисертації на здобуття наукового ступеня

доктора технічних наук

Київ - 2006

Дисертацією є рукопис.

Робота виконана в Національному технічному університеті України “Київський політехнічний інститут” Міністерства освіти і науки України на кафедрі автоматизації проектування енергетичних процесів і систем.

Науковий консультант - доктор технічних наук, професор

Верлань Анатолій Федорович,

Інститут проблем моделювання в енергетиці ім. Г.Є.Пухова НАН України, завідувач відділу моделювання динамічних систем

Офіційні опоненти: член-кореспондент НАН України, доктор технічних наук, професор метод моделювання високоградієнтний процес

Рєзцов Віктор Федорович,

Інститут відновлюваної енергетики НАН України, заступник директора з наукових питань;

доктор технічних наук, професор

Нестеренко Борис Борисович,

Інститут математики НАН України, заступник директора;

доктор технічних наук, професор

Бейко Іван Васильович,

Українсько-Угорський інститут кібернетики ім. Арпада Гьонца Міжрегіональної академії управління персоналом Міністерства освіти і науки України, директор.

Провідна установа Інститут кібернетики ім. В.М.Глушкова НАН України, відділ 265 - розподілених систем, що самоорганізуються.

Захист відбудеться 12 червня 2006 р. о 14-30 на засіданні спеціалізованої ради Д 26.002.02 у НТУУ “КПІ” (м. Київ, пр. Перемоги, 37, корп. 18, ауд. 306).

Відзиви на автореферат у двох екземплярах, завірені печаткою установи, просимо надсилати на адресу: 03056, м. Київ, пр. Перемоги, 37, вченому секретарю НТУУ “КПІ”.

З дисертацією можна ознайомитись в бібліотеці Національного технічного університету України “Київський політехнічний інститут”, м. Київ, пр. Перемоги, 37.

Автореферат розісланий “_10__”___05____ 2006 р.

Вчений секретар спеціалізованої ради, кандидат технічних наук, доцент М.М. Орлова

ЗАГАЛЬНА ХАРАКТЕРИСТИКА РОБОТИ

Актуальність теми. Математичне моделювання з використанням комп'ютерної техніки зараз є загальноприйнятим способом дослідження різноманітних природних явищ та створення нових технічних пристроїв. Комп'ютерні системи автоматизованого проектування скорочують терміни проектування та підвищують якість нових технічних систем. Значне місце в інженерній практиці займає моделювання об'єктів з розподіленими параметрами (ОРП), з якими стикаються при дослідженні процесів тепломасопереносу, руху рідини та газу, розповсюдження електромагнітних хвиль та інше. Теоретичною основою математичного моделювання ОРП є фундаментальні закони фізики: збереження енергії, маси та інші. Значну частину з них можна сформулювати математично у вигляді диференційних рівнянь у частинних похідних (ДРЧП), розв'язання яких за допомогою комп'ютерних систем (КС) вимагає великих витрат машинного часу, особливо у тривимірному випадку. Внаслідок необхідності вирішення все більш складних задач дослідження фізичних об'єктів, синтезу та проектування технічних систем, зокрема, проведення багатоваріантного їх аналізу з метою оптимізації параметрів, актуальною проблемою є скорочення машинних витрат та підвищення точності обчислень при комп'ютерній реалізації цих моделей. Основними методами розв'язання ДРЧП є методи скінчених різниць та скінчених елементів. Їх загальними рисами є заміна неперервної області визначення невідомої функції (температури, концентрації речовини та інше) дискретною множиною точок або елементів (сіткою), в яких і відшукується розв'язок. Відносно значень невідомої функції в вузлах цієї сітки формується та розв'язується алгебрична система рівнянь, що замінює диференційну задачу у кожний момент часу. Різноманітні варіанти цих методів відрізняються точністю, трудомісткістю, стійкістю та іншими характеристиками. Значний вклад у розвиток методів математичного моделювання ОРП та обчислювальних методів розв'язання ДРЧП внесли В.В.Ажогін, М.М.Бєляєв, А.Ф.Верлань, С.К.Годунов, В.С.Дейнека, В.Ф.Євдокімов, М.З.Згуровський, М.М.Каліткін, Л.О.Коздоба, В.І.Лаврик, І.І.Ляшко, В.Л.Макаров, Г.І.Марчук, Ю.М.Мацевитий, І.М.Молчанов, Б.Б.Нестеренко, Г.М.Положій, В.Л.Рвачов, В.Ф.Рєзцов, П.Роуч, О.А.Самарський, І.В.Сергієнко, В.В.Скопецький, О.О.Скоробагатько, А.М.Тихонов, Л.І.Турчак, С.Фарлоу, В.К.Хрущ та інші.

Етап дискретизації неперервної області визначення шуканої функції, тобто визначення кількості та розташування вузлів різницевої сітки, є дуже важливим і значною мірою визначає ефективність розв'язання задачі. При збільшенні кількості вузлів підвищується точність, але зростає розмірність алгебричних систем та, відповідно, час обчислень. Точність різницевого розв'язку залежить також від поведінки шуканої функції, а саме, зменшується у зонах великих її градієнтів. Отже, доцільно згущувати вузли у таких зонах та розташовувати їх рідко там, де функція змінюється досить плавно, що дозволяє скоротити розмірність систем алгебричних рівнянь без погіршення точності та час обчислень. Аналіз свідчить, що проблема створення адаптивних методів чисельного моделювання ОРП, орієнтованих на ефективну комп'ютерну реалізацію, є недостатньо вирішеною і потребує розвитку відповідних напрямків теорії обчислень, зокрема, напрямку формалізації та алгоритмізації способів побудови оптимальних різницевих сіток, що особливо актуально у випадках, коли зони великих градієнтів з часом змінюють своє положення у просторі (наприклад, рух променя лазера по поверхні матеріалу, розповсюдження забруднень в атмосфері або водоймищі).

Важливою проблемою, що виникає при побудові різницевої сітки, є криволінійність границь області зміни аргументів, що ускладнює алгоритм і досить часто зустрічається на практиці. Методам побудови сіток, що адаптовані до криволінійних границь та поведінки розв'язку, присвячені праці Б.М.Азаренка, Д.Андерсона, П.М.Вабищевича, М.О.Дар'їна, Л.М.Дегтярьова, В.В.Дроздова, Х.А.Дуайера, С.А.Іваненка, В.Д.Лисейкіна, В.І.Мажукіна, Г.П.Прокопова, Дж.Таннехілла, Д.Ф.Томпсона, К.Флетчера та інших. Істотним недоліком існуючих методів є недостатня формалізація та автоматизація етапу підготовки задачі до розв'язання, що полягає, наприклад, у необхідності ручного вибору вагових функцій, визначення початкового положення граничних вузлів та інше. Існуючі адаптивні методи призначені, головним чином, для зменшення похибки обчислень, а не для забезпечення можливостей одержання розв'язку з наперед заданою точністю.

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

Зв'язок з науковими програмами, планами, темами. Дисертаційна робота виконувалася відповідно до пріоритетних напрямків наукових досліджень, державних програм і тем, зокрема:

- науково-технічної програми “Розробка науково-методичних основ системи прогнозування генетичного ризику впровадження нових технологій та забруднення навколишнього середовища ГРАНІТ” (номер державної реєстрації теми УКРІНТЕІ №080999 Р0005), яка розроблена на виконання Указу Президента України від 17 січня 1995 року № 53/95 "Про систему прогнозування генетичного ризику впровадження нових технологій та забруднення навколишнього середовища". Запропонований автором адаптивний метод був застосований для моделювання розповсюдження забруднень в атмосфері;

- теми "Структурно-алгоритмічні методи й засоби комп'ютерного моделювання складних енергетичних об'єктів з розподіленими й змінними параметрами" (номер державної реєстрації теми 0101U000024). Автором розроблено метод моделювання процесів перенесення тепла у теплових трубах;

- теми “Розробка інформаційної технології моделювання і розрахунку теплових режимів перспективних конструкцій теплонавантажених елементів та пристроїв засобів обчислювальної техніки та керування” (номер державної реєстрації теми 0104U003807). Запропонований автором адаптивний метод був використаний для визначення температурних полів теплопередаючих пристроїв;

- теми “Створення спеціальних низькотемпературних теплових труб як складової комбінованої системи охолодження інфрачервоної техніки виробів мобільного призначення” (номер державної реєстрації теми 0202U004040). Автором розроблений адаптивний алгоритм для моделювання температурних полів контурних теплових труб.

Мета і завдання дослідження. Мета роботи полягає у створенні, дослідженні і застосуванні адаптивних методів чисельного моделювання високоградієнтних процесів в об'єктах з розподіленими параметрами з орієнтацією на ефективну комп'ютерну реалізацію і забезпеченням при цьому автоматизованої побудови оптимальних різницевих сіток, одержання результатів моделювання з мінімальними витратами машинних ресурсів та апріорі заданою точністю. Відповідно до зазначеної мети поставлені такі задачі:

Проаналізувати відомі підходи та методи моделювання високоградієнтних процесів в ОРП та виявити фактори, які істотно впливають на його ефективність.

Розробити метод оцінки похибки результату, який дозволить одночасно обчислити її чутливість до змін кроків сітки і скоротить додаткові розрахунки на визначення характеру трансформації сітки .

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

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

Розробити методику урахування нерівномірності сітки при формуванні різницевих рівнянь для одно-, дво- та тривимірних задач, що дасть змогу скоротити кількість вузлів у зонах плавної зміни шуканої функції.

Розробити алгоритмічне та програмне забезпечення для моделювання високоградієнтних процесів в ОРП з використанням адаптивних методів та дослідити їх ефективність при розв'язанні тестових і прикладних задач різної розмірності.

Виявити особливості використання адаптивних методів для моделювання процесів з вузьколокалізованими зонами зовнішнього впливу великої потужності та розробити алгоритми побудови різницевої сітки для таких випадків.

Сформулювати рекомендації щодо використання різних типів адаптивних сіток в залежності від властивостей процесу, що моделюється.

Об'єкт дослідження - математичні моделі високоградієнтних процесів в ОРП.

Предмет дослідження - адаптивні чисельні методи моделювання ОРП, орієнтовані на ефективну реалізацію комп'ютерними засобами.

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

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

1. Запропоновано й обґрунтовано класифікацію нестаціонарних моделей ОРП з погляду наявності зон з великими градієнтами й динаміки їхніх координат і сформульовані рекомендації по застосуванню адаптивних сіток різних типів в залежності від класу моделей, що скорочує час обчислень.

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

3. Запропоновано й обґрунтовано метод оцінки похибки розв'язку багатовимірних ДРЧП, який заснований на використанні вкладених сіток і на відміну від відомих методів, одночасно з похибкою визначає її чутливість до зміни конфігурації різницевої сітки, що скорочує додаткові обчислення на визначення характеру трансформації сітки.

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

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

6. Для розв'язання нежорстких моделей ОРП вперше запропоновано метод побудови змінної рівномірної різницевої сітки, що дозволяє одержувати результат з наперед заданою точністю і скорочує витрати на побудову адаптивної сітки.

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

8. Удосконалено метод розв'язання систем різницевих рівнянь, які формуються у методах розщеплення, що скорочує необхідний об'єм пам'яті КС.

9. Для урахування нерівномірності сітки в розрахунковій області удосконалено методи формування різницевих рівнянь.

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

Практичне значення одержаних результатів визначається тим, що запропоновані методи чисельного моделювання високоградієнтних процесів в ОРП дозволяють значно розширити клас математичних моделей щодо ефективного комп'ютерного дослідження багатьох видів фізичних процесів та створення сучасних зразків об'єктів нової техніки, значно скорочувати терміни розрахунків і досліджень, підвищувати та контролювати точність результатів моделювання процесів в ОРП у різних галузях науки і техніки. Теоретичні і практичні результати дисертації були використані в наукових дослідженнях при виконанні низки бюджетних і договірних тем, зокрема, при моделюванні процесів поширення забруднювачів в атмосфері, що виконувалося в НТУУ “КПІ” в рамках НТП ГРАНІТ (акт про впровадження від 9.03.06); при моделюванні впливу лазерного випромінювання на матеріали, що проводилося кафедрою лазерних технологій НТУУ “КПІ” (довідка про впровадження від 9.03.06); при моделюванні температурних полів контурних теплових труб, що виконувалося інженерним центром “Перспектива” НТУУ “КПІ” (довідка про впровадження від 17.01.05). Також результати дисертації використовуються в навчальному процесі НТУУ “КПІ” при викладанні курсів лекцій “Чисельні методи в інформатиці”, “Чисельні методи в задачах математичної фізики”, у курсовому і дипломному проектуванні студентів спеціальностей “Інформаційні технології проектування”, “Програмне забезпечення автоматизованих систем” (довідка про впровадження від 9.03.06) і відображені у навчальних посібниках “Математичне забезпечення САПР” (Лук'яненко С.О. - К. : ІСДО, 1996. -148с.), “Основи обчислювальних методів розв'язування диференціальних рівнянь” (Лук'яненко С.О. - К. : ІСДО, 1998. -212с.) із грифом Міністерства освіти і науки України.

Особистий внесок здобувача в опублікованих працях полягає в обґрунтуванні авторського підходу до побудови адаптивних методів моделювання ОРП і розробці адаптивних методів, які забезпечують наперед задану точність та мінімізують витрати машинного часу за рахунок скорочення кількості вузлів та їх оптимального розташування. Особисто автором розроблені обчислювальні модулі програмного забезпечення для розв'язання багатовимірних задач моделювання ОРП, що використовує запропоновані адаптивні методи.

Сформульовані в дисертації наукові результати, висновки і рекомендації належать особисто автору і є його науковим доробком.

Апробація результатів дисертації. Основні положення дисертації доповідались, обговорювались і були схвалені на постійно діючому науковому семінарі з міжнародною участю “Економічні проблеми енергетичної безпеки” за темою ”Еколого-економічні аспекти інноваційного розвитку паливно-енергетичного комплексу” (м.Трускавець, 27-28 листопада 2003р.), науково-практичній конференції “Оцінка ризику захворювань від забруднення навколишнього середовища та впровадження нових технологій” (м. Київ, 24 грудня 2003р.), другій міжнародній конференції “Обчислювальна та прикладна математика” (м. Київ, 27 вересня 2004р.), науковому семінарі з міжнародною участю “Інформаційні технології в забезпеченні економічної безпеки держави” (м. Біла Церква, 27-29 квітня 2005р.), сьомому міжнародному семінарі “Сучасні інформаційні технології” (м. Мінськ, 2-9 липня 2005р.), міжнародній конференції “Інформаційні технології в управлінні енергетичними системами” (м. Київ, 18-19 жовтня 2005р.).

Публікації. Наукові положення, висновки і рекомендації дисертаційного дослідження опубліковані в 30 роботах, з яких 27 відповідають вимогам ВАК України, у тому числі: в 1 одноосібній монографії, 11 - у монографічних виданнях, 15 - у наукових журналах та збірниках наукових праць; 3 - публікації матеріалів конференцій.

Структура й обсяг дисертації. Дисертація складається зі вступу, шести розділів, висновків, списку використаних джерел з 210 найменувань, додатків і містить 50 таблиць та 98 ілюстрацій. Повний обсяг дисертації складає 315 сторінок, з яких список використаних джерел займає 20 сторінок, додатки - 6 сторінок.

ОСНОВНИЙ ЗМІСТ РОБОТИ

У першому розділі “Проблеми та методи математичного моделювання об'єктів з розподіленими параметрами на основі скінчено-різницевого підходу” розглядається класифікація моделей об'єктів та виділяється їх важливий клас: моделі мікрорівня проектування, або моделі об'єктів з розподіленими параметрами (МОРП), що у багатьох випадках мають вигляд диференційних рівнянь у частинних похідних (ДРЧП). Вони широко використовуються не тільки при проектуванні нових технічних систем, але і при дослідженні різноманітних природних явищ у теплофізиці, екології, аерогідродинаміці та ін. Невідомими функціями при цьому є, наприклад, температура, тиск, концентрація та інші характеристики просторово-часових фізичних полів, а їх аргументами - просторові координати та час. Аналітичні методи розв'язання таких задач існують лише для деяких їх класів, тому основою сучасного підходу до розв'язання МОРП виступають чисельні методи з використанням обчислювальної техніки. Але ці методи є дуже трудомісткими, особливо для тривимірних задач, що ускладнює їх використання в ітераційних процедурах параметричного синтезу об'єктів нової техніки.

На складність розв'язання МОРП, крім інших чинників, значно впливає вигляд розв'язку, а саме наявність зон з великими градієнтами, що погіршує точність апроксимації диференційної задачі. З цієї точки зору можна провести класифікацію МОРП за критеріями наявності великих градієнтів, затухання їх з часом та зміни положення у просторі. Назвемо МОРП жорсткою, якщо у її розв'язку присутня хоча б одна зона з великими градієнтами. Якщо великі градієнти розв'язку присутні лише у деякому діапазоні модельного часу (затухання або виникнення неоднорідності), а в інші моменти часу розв'язок є гладким, МОРП назвемо частково жорсткою, інакше - повністю жорсткою. Якщо зона великих градієнтів змінює з часом своє положення у просторі, МОРП назвемо жорсткою динамічною (наприклад, модель руху променя лазера по поверхні матеріалу), а інакше - статичною.

Існує досить багато обчислювальних методів розв'язання ДРЧП, які різняться способом заміни диференційної задачі набором алгебричних. Умовно їх розділяють на дві великі групи методів: скінчених різниць та скінчених елементів. Їх основними характеристиками є збіжність, точність апроксимації, стійкість та інше. Незважаючи на їх різноманіття, вони мають деякі спільні риси. При їх використанні неперервну область зміни аргументів шуканої функції замінюють кінцевою (дискретною) множиною точок (вузлів), які визначають різницеву сітку; на сітці замість функцій неперервного аргументу розглядають функції дискретного аргументу, визначені в її вузлах, або використовують їх локальну апроксимацію за допомогою функцій, визначених на обмежених вузлами елементах; диференційне рівняння для кожного моменту модельного часу замінюють системою алгебричних рівнянь (різницевими рівняннями), початкові й крайові умови також замінюють різницевими початковими й крайовими умовами.

У результаті одержують різницеву крайову задачу у вигляді однієї або декількох систем алгебричних рівнянь. Розв'язком різницевої задачі є сіткова функція, визначена у вузлах сітки, тобто на дискретній множині точок, або кусково-неперервна функція для методу скінчених елементів.

Складність та трудомісткість розв'язання багатовимірних задач моделювання ОРП обумовлюють актуальність пошуку шляхів підвищення ефективності їх чисельного розв'язання, причому способом досягнення цієї мети у даній роботі прийнято скорочення кількості рівнянь у алгебричному аналогу диференційного рівняння за рахунок раціонального розташування вузлів сітки.

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

В залежності від форми чарунок сітки діляться на трикутні, чотирикутні і т.д. В залежності від типу ліній або поверхонь, що обмежують чарунки, сітки можуть бути криволінійними або прямолінійними. Якщо розміри усіх чарунок однакові, сітки називають рівномірними, інакше - нерівномірними. Якщо з часом розміри чарунок не змінюються, сітки називають фіксованими, інакше - змінними. Основними проблемами при дискретизації області є урахування криволінійності границь та адаптація сітки до поведінки функції. У багатьох практичних випадках границя області має криволінійну форму і граничні точки області не збігаються з вузлами сітки. Тоді при формулюванні різницевої задачі можна: ввести додаткові вузли на перетині ліній x=const та y=const з границею; границю приблизно замінити ламаною, яка проходить через близькі до границі вузли і на цю ламану перенести граничні умови; виконати перетворення координат, при якому фізична область неправильної форми трансформується у розрахункову область у вигляді прямокутника або паралелепіпеда відповідно для дво- або тривимірних задач.

В останньому випадку задача полягає в знаходженні відображення , що переводить вузли сітки з фізичної області G у розрахункову область . При цьому вимагають, щоб відображення було взаємно однозначним; лінії сітки у фізичній області повинні бути гладкими; чарунки сітки у фізичній області не повинні бути занадто перекошені, оскільки це збільшує похибку; сітка у розрахунковій області повинна бути рівномірною. Методи перетворення координат при побудові розрахункових сіток поділяють на такі групи: методи теорії функцій комплексної змінної (можуть бути застосовані лише для двовимірних задач); алгебричні методи; диференційні методи, що базуються на розв'язуванні допоміжних ДРЧП; варіаційні методи.

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

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

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

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

вибрати нові величини кроків по часу та просторовим координатам так, щоб при максимальній величині кроків похибка знаходилась в заданих межах.

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

Отже, основними етапами запропонованих адаптивних методів на кожному часовому кроці є:

обчислення розв'язку на декількох вкладених сітках;

використання одержаної інформації для оцінки похибки та визначення характеру трансформації сітки для кожного вузла;

побудова нової сітки на основі її функцій трансформації;

апроксимація розв'язку у вузлах нової сітки.

Нехай для розв'язання диференційної задачі використано деякий різницевий метод з точністю апроксимації на рівномірних сітках

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

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

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

Нехай розв'язок для моменту часу вже знайдений i треба знайти його для моменту часу , причому перевiрити його точнiсть та вибрати нову сiтку для переходу на наступний шар . Сiтку для шару вважаємо вiдомою. Якщо за допомогою рiзницевої схеми перейти з шару на шар з базовим кроком та просторовими кроками, що задаються вiдомою нерiвномiрною сiткою, то для вузла з координатами виконується рiвнiсть (2). Якщо повторно перейти з шару k на шар k+1 з кроком , двiчi застосувавши вiдповiдну рiзницеву схему, то одержимо результат , для якого виконується рiвнiсть

Наступнi n рiвнянь шуканої системи одержимо, якщо повторимо перехiд з k-го на (k+1)-й шар почергово з удвiчi щiльнiшою сiткою спочатку в напрямку x₁, потiм x₂ i так далi. Одержанi в кожному випадку для кожного вузла значення функцiї позначимо через , і т.д. В результатi для кожного вузла сiтки на (k+1)-му шарi одержимо систему лiнiйних алгебричних рiвнянь розмірністю вiдносно :

Із цієї системи для кожного вузла визначаються коефіцієнти . У тривимірному випадку:

Це дозволяє для кожного вузла за формулою (2) визначити похибку:

Максимальне її значення по усім вузлам є локальною похибкою на (k+1)-му шарi і вона порівнюється з допустимою похибкою . Якщо , то результати анулюються і на базi k-го шару знову шукаємо розв'язок на ущільненій сітці. Якщо , крок приймається і уточнюється розв'язок:

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

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

Для цього необхідно у кожному вузлі розв'язати рівняння

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

що дає систему з (n+1)-го рівняння з (n+1)-ю невідомою величиною. Кожен набір з коефіцієнтів діє в околі свого вузла, кожен коефіцієнт дає допустиму зміну середнього кроку у даному координатному напрямку.

Спочатку з першого рівняння системи (4) знайдемо для кожного вузла коефіцієнт зміни часового кроку : .

Оскільки обчислення значення невідомої функції ведеться для одного й того ж моменту часу t_k+2, цей коефіцієнт повинен бути однаковим для усіх вузлів даного часового шару, тому з величин виберемо мінімальну по усім вузлам: _. При цьому допустима похибка, що припадає на інші доданки формули (3), збільшиться і її знайдемо за формулою

Позначимо праву частину цієї рівності через e₁ і знову вважатимемо, що складові похибки по координатам повинні бути однаковими:

З першого рівняння цієї системи знайдемо для кожного вузла сітки коефіцієнт зміни кроку h₁ по першій просторовій змінній:. Для кожного вузла першої просторової координати знайдемо функцію трансформації сітки по цій координаті:. Побудуємо нерівномірну сітку у напрямку х₁ за допомогою алгоритму, викладеному далі. Після цього стануть відомими відкориговані значення коефіцієнтів в_1,i , які не перевищують раніше обчислених. За рахунок цього збільшиться допустима похибка, що припадає на доданки формули (5), починаючи з другого. Їх суму

позначимо через e₂ і знову вимагатимемо рівності складових похибки:

Обчислимо для кожного вузла , знайдемо функцію трансформації та побудуємо нерівномірну сітку у напрямку х₂ і т.д. Після повної побудови нової нерівномірної сітки на (k+1)-му шарі засобами багатовимірної інтерполяції визначаються значення у нових вузлах цього опорного шару і починається спочатку перехід на наступний (k+2)-й шар.

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

Оскільки сума кроків нової сітки повинна дорівнювати довжині відрізку зміни координати і, крім того, на кроки накладаються додаткові умови, не можна безпосередньо використати для формування сітки. Побудова сітки по одній з координат полягає в знаходженні мінімуму функціонала при умовах:, ,.

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

Якщо взяти х_1,нов малим, то попадемо в область, де допускається значно більший крок і одержаний результат буде дуже точним, але суттєво збільшаться витрати машинного часу. Таким чином, необхідно вибрати максимальний х_1,нов такий, щоб його значення було меншим або рівним допустимого кроку h_доп(х) на відрізку [x₀,x_1,нов]:. Наступний вузол х_2,нов знаходиться аналогічно:. У загальному випадку:.

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

У третьому розділі “Адаптивні методи вкладених сіток для одновимірних моделей нестаціонарних процесів” розроблено алгоритми моделювання одновимірних процесів за допомогою двох запропонованих адаптивних методів: з використанням змінних рівномірних (ЗРС) та змінних нерівномірних (ЗНС) сіток - та досліджено їх ефективність. При використанні ЗРС значно спрощується алгоритм: після визначення коефіцієнтів трансформації сітки у деякому координатному напрямку не потрібно визначати функцію трансформації , функцію допустимого кроку та будувати нерівномірну сітку у даному напрямку. Замість цього необхідно знайти коефіцієнт та обчислити величину нового кроку для даної координати. Це зменшує ті додаткові витрати машинного часу, які необхідні, крім формування та розв'язання різницевих рівнянь, для реалізації властивості адаптивності. Але такий підхід доцільно використовувати лише для нежорстких моделей, оскільки він не забезпечує скупчення вузлів у зоні великих градієнтів у випадку їх наявності.

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

У початковий момент часу температура стержня дорівнює одиниці і потім його кінці миттєво охолоджуються до температури 0С. Математична постановка:

Задача має розривні крайові умови та великі значення похідних на кінцях відрізку [0,1], особливо на початку досліджуваного часового інтервалу.

Якщо розв'язувати цю задачу з фіксованою рівномірною сіткою (ФРС), це потребує розрахунку з досить малими кроками по обом змінним та великих витрат машинного часу, тобто поведінка функції U(x,t) для малих значень t у цьому випадку вимагає густої сітки для усієї розрахункової області.

На протилежність цьому, адаптивний алгоритм зі ЗРС виконує обчислення з малими кроками лише в початковий період, а далі істотно збільшує їх значення. Оскільки з ростом часу градієнт функції для цієї задачі зменшується, то зростає і крок по часу. Але тут при фіксованому t вузли розташовані рівномірно і не враховується плавна зміна функції у середній частині відрізку [0,1], де вузли можна розташувати значно рідше. Великий градієнт на кінцях цього відрізку приводить при ЗРС до малого кроку по х на усьому відрізку. Цього недоліку позбавлена ЗНС (рис.5в), для якої просторові вузли скупчуються лише на кінцях відрізку, що скорочує розмір системи алгебричних рівнянь та витрати машинного часу. Порівняння ефективності цих трьох видів сіток показує, що адаптивний алгоритм зі ЗНС скорочує час розрахунку порівняно з використанням алгоритму, що базується на фіксованій сітці, у 3-10 разів і порівняно з використанням ЗРС - у 3-4 рази.

Порівняння часу розрахунків одновимірної тестової задачі №1 з фіксованою рівномірною, змінною рівномірною та змінною нерівномірною сітками

У другій тестовій задачі моделюється розповсюдження одновимірної теплової хвилі. Цей процес описується квазілінійним рівнянням теплопровідності:

При порівнянні процесу розв'язання цієї задачі з трьома типами сіток одержані показники ефективності аналогічні попередньому.

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

Проведені дослідження показують, що ЗРС, у яких просторовий крок сітки однаковий для кожного моменту часу і змінюється при переході на наступний часовий шар, а також є змінним часовий крок, є проміжним засобом підвищення ефективності моделювання об'єктів з розподіленими параметрами. Такі сітки вже дають змогу одержувати результат з наперед заданою точністю. Для задач з гладким розв'язком (нежорстких задач) можна підібрати фіксовану рівномірну сітку (ФРС), яка забезпечить достатню точність з малими витратами машинного часу і адаптивні сітки не дають переваг, окрім можливості автоматичного досягнення наперед заданої точності. Для задач, розв'язок яких має зони з великими градієнтами, що не змінюють з часом своє положення у просторі та величину (статичні жорсткі задачі), можна підібрати фіксовану нерівномірну сітку (ФНС), яка теж забезпечить достатню точність з малими витратами машинного часу і адаптивні сітки не дають переваг, окрім можливості автоматичного досягнення наперед заданої точності і позбавлення необхідності вручну будувати нерівномірну сітку. Для частково жорстких задач, розв'язок яких має зони з великими градієнтами, що не змінюють з часом своє положення у просторі, але величина градієнтів істотно змінюється (гладкий розв'язок перетворюється у негладкий або навпаки) ФНС є неефективною за витратами машинного часу і доцільно використовувати найпростішу форму адаптивних сіток - ЗРС. Для динамічних жорстких задач, розв'язок яких має зони з великими градієнтами, що змінюють з часом своє положення у просторі (моделювання руху теплової хвилі або променя лазера), неефективними за витратами машинного часу є не тільки фіксовані сітки, але й ЗРС, тому необхідно використовувати ЗНС. Таким чином, використання адаптивних методів ускладнює програму та вимагає додаткових розрахунків на оцінку поведінки розв'язку, перебудову сітки, але скорочує розмірність систем алгебричних рівнянь. Адаптивний алгоритм зі ЗНС скорочує час розв'язання одновимірних задач порівняно з використанням алгоритму, що базується на фіксованій сітці, у 3-10 разів і порівняно з використанням ЗРС - у 2-4 рази.

Четвертий розділ “Чисельне моделювання двовимірних нестаціонарних процесів” присвячений розробці, обґрунтуванню та дослідженню алгоритмів розв'язання адаптивними методами лінійних та нелінійних двовимірних задач. Для дослідження ефективності методів порівнюються результати розрахунків тестових задач методами змінних напрямків та розщеплення по координатах з використанням ЗРС та ЗНС.

Перша тестова задача моделює зміну температурного поля пластини при охолодженні двох її сторін. Задача має простий точний розв'язок і тому зручна для дослідження ефективності методу. Математична постановка:

Її точним розвязком є функція U(x,y,t)=e^-2xyt , яка у початковий момент часу являє собою горизонтально розташовану площину U=1.

Для розв'язання даної задачі застосовано метод змінних напрямків, модифікований для нерівномірної сітки. Часовий крок виконується у два етапи:

При обчисленнях відповідно до збільшення градієнтів розв'язку алгоритм зменшує крок за часом та просторові кроки (рис.8).

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

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

Порівняння двох типів адаптивних сіток показало, що для ЗРС скупчення вузлів відбувається по усій області, а для ЗНС - лише на початку координат, що зменшує загальну кількість вузлів і, відповідно, час розрахунку.

Порівняння ефективності розрахунку двовимірної тестової задачі №1 методом змінних напрямків з фіксованою та змінною нерівномірною сітками

Дослідження ефективності розв'язання даної задачі методом змінних напрямків та методом розщеплення по координатах зі змінною нерівномірною сіткою показало перевагу методу змінних напрямків.

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

ільки градієнти функції з ростом t зменшуються, алгоритм збільшує величини кроків сітки.

Порівняння алгоритмів з фіксованими та змінними кроками показує, що для цієї задачі адаптивний метод зі змінною нерівномірною сіткою скорочує витрати машинного часу приблизно у 10 разів.

Порівняння часу розрахунку двовимірної тестової задачі №2 методом змінних напрямків з фіксованою та змінною нерівномірною сітками

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

Таким чином, застосування адаптивних методів вкладених сіток для лінійних та нелінійних двовимірних задач забезпечує скорочення витрат машинного часу у 5-10 разів порівняно з фіксованою рівномірною сіткою при однаковій точності.

У п'ятому розділі “Методи та алгоритми моделювання тривимірних процесів” розроблено та досліджено алгоритми адаптивних методів для розв'язання тривимірних задач на нерівномірній сітці двоциклічним шостиетапним методом розщеплення по координатах: де ₁, ₂, ₃ - різницеві оператори:

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

Експериментально за допомогою тестової задачі проведено дослідження ефективності адаптивних методів, яке показало, що при однаковій точності результату вони працюють у 4-5 разів швидше.

У шостому розділі “Комп'ютерні засоби та розв'язання прикладних задач” розглянуто застосування адаптивних методів вкладених сіток до розв'язання складних прикладних задач моделювання розповсюдження забруднювачів в атмосфері, переносу тепла у теплопередаючих пристроях, впливу лазерного випромінювання на матеріали, а також структуру програмного комплексу моделювання ОРП.

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

У першій задачі для моделювання розповсюдження забруднювачів на велику відстань (більше 100км) використано осереднене по висоті значення концентрації забруднювачів U(x,y,t). Нехай внаслідок вибуху на деякій висоті утворилася сферична хмара, заповнена газоаерозольною речовиною з відомою концентрацією. Необхідно дослідити зміну концентрації речовини з часом на цій висоті внаслідок дії вітрового потоку та механізму дифузії. Задача зводиться до розв'язання двовимірного рівняння нестаціонарної турбулентної дифузії

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

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

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

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

Порівняння часу розрахунку тестової задачі з фіксованою та змінною нерівномірною сітками

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