Проанализированный документ: Kozub_G.pdf Лицензия: ВОЛОДИМИР МАТІЄВСЬКИЙ
Детальный анализ тела документа:
Диаграмма соотношения частей:
Детали обработанных ресурсов:
162 - ОК /
5 - Ошибок
Активные ссылки (URL-адреса, извлеченные из документа):
Детальный анализ документа:
МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ ДЕРЖАВНИЙ ЗАКЛАД
«ЛУГАНСЬКИЙ НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ ІМЕНІ ТАРАСА ШЕВЧЕНКА»
Навчально-науковий інститут фізики, математики та інформаційних технологій Кафедра фізико-технічних систем та інформатики Козуб Галина Олександрівна ЧИСЕЛЬНЕ МОДЕЛЮВАННЯ ФІЗИКО-МЕХАНІЧНИХ ПОЛІВ В ТВЕРДИХ ТІЛАХ кваліфікаційна робота здобувача вищої освіти другого (магістерського) рівня освітньої програми
«Комп’ютерні науки та інформаційні технології»
за спеціальністю 122 Особистий підпис __________ Науковий керівник _____________ Ольга СМАГІНА, кандидат педагогічних наук, доцент кафедри фізико-технічних систем та інформатики Завідувач кафедри ______________ Юрій КОЗУБ, доктор технічних наук, доцент кафедри фізико-технічних систем та інформатики Полтава – 2023 АНОТАЦІЯ Козуб Г.О. Тема: Чисельне моделювання фізико-механічних полів в твердих тілах Спеціальність: 122 Установа: ДЗ ЛНУ імені Тараса Шевченка, 2023 р. Кваліфікаційна робота містить: 90 стор., 4 табл., 15 рис., 43 джерело, 2 додатки. Об’єкт дослідження – процеси термопружного деформування елементів конструкцій з твердих тіл. Предмет дослідження – величини параметрів напружено-деформованого стану конструкцій в умовах термосилового навантаження. Мета роботи – розв’язання проблеми визначення термомеханічних полів, напружень та деформацій в просторових конструкціях із анізотропних матеріалів. Методи дослідження: Теоретичні методи: аналіз науково-технічних джерел з проблем дослідження. Емпіричні методи: аналіз проблеми чисельного моделювання фізико-механічних полів в твердих тілах; дискретизація твердих тіл чисельними методами. Експериментальні методи: тестування розробленого додатку на експериментальних задачах. Результати роботи. Досліджено основні вихідні і розрахункові співвідношення методу скінченних елементів для задач теплопровідності конструкцій, а також алгоритми їх розв’язання. Проведено чисельні дослідження для обґрунтування достовірності результатів. Розроблено програмний комплекс для чисельного розв’язання зв’язаних задач термопружності з використанням МСЕ. Ключові слова: МСЕ, НДС, ТЕРМОПРУЖНІСТЬ, ТЕПЛОПРОВІДНІСТЬ, МІРЕЛА+ АBSTRАСT Kоzub H.О. Thеmе: Numеrісаl mоdеlіng оf рhуsісаl аnd mесhаnісаl fіеlds іn sоlіd bоdіеs Sресіаltу: 122 Іnstіtutіоn: Luhаnsk Tаrаs Shеvсhеnkо Nаtіоnаl Unіvеrsіtу, 2023 Quаlіfісаtіоn wоrk соntаіns: 90 раgеs, 4 tаblеs, 15 fіgurеs, 43 sоurсеs, 2 арреndісеs. Оbjесt оf rеsеаrсh рrосеssеs оf thеrmоsрrіng dеfоrmаtіоn оf struсturаl еlеmеnts mаdе оf sоlіd bоdіеs. Subjесt оf rеsеаrсh - vаluеs оf thе раrаmеtеrs оf thе strеss-strаіn stаtе оf struсturеs undеr соndіtіоns оf thеrmоfоrсе lоаdіng Рurроsе оf thе studу sоlvіng thе рrоblеm оf dеtеrmіnіng thеrmоmесhаnісаl fіеlds, strеssеs аnd dеfоrmаtіоns іn sраtіаl struсturеs mаdе оf аnіsоtrоріс mаtеrіаls. Rеsеаrсh mеthоds: Thеоrеtісаl mеthоds: аnаlуsіs оf sсіеntіfіс аnd tесhnісаl sоurсеs оn rеsеаrсh рrоblеms. Еmріrісаl mеthоds: аnаlуsіs оf thе рrоblеm оf numеrісаl mоdеlіng оf рhуsісаl аnd mесhаnісаl fіеlds іn sоlіd bоdіеs; dіsсrеtіzаtіоn оf sоlіd bоdіеs bу numеrісаl mеthоds. Ехреrіmеntаl mеthоds: tеstіng thе dеvеlореd аррlісаtіоn оn ехреrіmеntаl tаsks. Rеsults оf rеsеаrсh. Thе mаіn іnіtіаl аnd саlсulаtіоn rаtіоs оf thе fіnіtе еlеmеnt mеthоd fоr рrоblеms оf thеrmаl соnduсtіvіtу оf struсturеs, аs wеll аs аlgоrіthms fоr thеіr sоlutіоn, wеrе studіеd. Numеrісаl studіеs wеrе соnduсtеd tо substаntіаtе thе rеlіаbіlіtу оf thе rеsults. А sоftwаrе соmрlех hаs bееn dеvеlореd fоr thе numеrісаl sоlutіоn оf соuрlеd рrоblеms оf thеrmоеlаstісіtу usіng FЕM. Kеуwоrds: FЕM, SDS, THЕRMОЕLАSTІСІTУ, THЕRMАL СОNDUСTІVІTУ, MІRЕLА+ ЗМІСТ Стор. ВСТУП ..................................................................................................................... 6 РОЗДІЛ 1. АНАЛІЗ СУЧАСНОГО СТАНУ ІСНУЮЧИХ РІШЕНЬ ............... 10 1.1. Проблеми та задачі термомеханіки твердих тіл ......................................... 10 1.2. Сучасні підходи до розв’язання задач термомеханіки у тривимірній постановці .............................................................................................................. 18 1.3. Висновки до розділу 1 ................................................................................... 26 РОЗДІЛ 2. МЕТОДИКА РОЗВ’ЯЗАННЯ ЗАДАЧ ТЕПЛОПРОВІДНОСТІ КОНСТРУКЦІЙ З ТВЕРДИХ ТІЛ ....................................................................... 28 2.1. Скінченний елемент для розв’язання задач теплопровідності .................. 28 2.2. Нестаціонарні задачі теплопровідності композитів ................................... 40 2.3. Дослідження збіжності рішення ................................................................... 46 2.4. Висновки до розділу 2 ................................................................................... 52 РОЗДІЛ 3. ПРОГРАМНА РЕАЛІЗАЦІЯ ТА МЕТОДИКА РОЗВ’ЯЗАННЯ ЗАДАЧ ТЕРМОПРУЖНОСТІ ............................................................................. 55 3.1. Структура обчислювального комплексу ............................... 55 3.2. Підсистеми розв’язування задач теплопровідності ..................... 57 3.3. Формування вихідних даних ......................................................................... 60 3.4. Програмне забезпечення підсистеми ............................................. 61 3.5. Побудова алгоритму розв'язання задачі термопружності ......................... 66 3.6. Програмне забезпечення підсистеми ..................................... 70 3.7. Висновки до розділу ...................................................................................... 72 Список використаних джерел .............................................................................. 75 ДОДАТКИ .............................................................................................................. 80 Додаток А. Фрагмент коду програми модулю ................................. 80 Додаток Б. Сертифікат впровадження ................................................................ 90 ПЕРЕЛІК УМОВНИХ ПОЗНАЧЕНЬ, СИМВОЛІВ, СКОРОЧЕНЬ МСЕ – метод скінчених елементів; МСР – метод скінченних різниць; ВРМ – варіаційно-різницевий метод; ММК – метод Монте-Карло; МП – метод прямих; МГЕ – метод граничних елементів; МІРЕЛА – міцність і руйнування еластомірних матеріалів; БД – база даних; НДС – напружено-деформований стан; СЕ – скінчений елемент; САD – Соmрutеr Аіdеd Dеsіgn; САЕ – Соmрutеr Аіdеd Еngіnееrіng; RАD – Rаріd Оf Арреndіхеs Dеvеlорmеnt; ПЗ – програмне забезпечення; ІDЕ – Іntеgrаtеd Dеsіgn Еnvіrоnmеn. 6 ВСТУП Актуальність теми. Сучасні досягнення в галузі машинобудування, при проектуванні будівельних конструкцій, аерокосмічної і військової техніки тісно пов’язані з використанням еластомерів як конструкційних матеріалів. Елементи конструкцій з таких матеріалів, як правило, працюють в умовах інтенсивного термосилового навантаження та супроводжуються ефектами як зворотної (термопружні), так і незворотної дії (розігрів внаслідок внутрішньої дисипації механічної енергії). Високі вимоги до надійності і економічності конструкцій, їх безвідмовної роботи ведуть до необхідності розгляду складних просторових задач термопружного деформування при наявності різноманітних супутніх факторів, що ускладнюють аналіз їх поведінки. Серед таких факторів можна виділити набуті або запроектовані початкові напруження, анізотропію властивостей матеріалів, які значною мірою впливають на експлуатаційні характеристики конструкцій, що розглядаються. Опис поведінки матеріалу з урахуванням теплових ефектів проводиться в рамках термомеханіки. Розв’язання таких задач можливе лише на основі чисельних методів, серед яких провідне місце займає метод скінченних елементів (МСЕ), що засвідчено його широким використанням в інженерній практиці. Таким чином, розробка та реалізація об’єктно-орієнтованих ефективних підходів до розв’язання просторових задач термопружного деформування конструкцій в механіки твердих тіл, їх термомеханічної поведінки є актуальною науковою проблемою сучасної механіки та представляє практичний інтерес. Мета роботи є розв’язання проблеми визначення термомеханічних полів, напружень та деформацій в просторових конструкціях із анізотропних матеріалів. Поставлена мета досягається вирішенням наступних завдань: 7 - створення чисельної методики моделювання термопружного деформування попередньо навантажених конструкцій; - розробка модифікованих розрахункових співвідношень просторових скінченних елементів для розв’язання задач теплопровідності анізотропних тіл; - побудова алгоритмів чисельного моделювання термомеханічних полів в твердих тілах на основі розв’язання задачі термопружності; - створення пакету прикладних програм
для визначення напружено- деформованого стану конструкцій із композитів; - аналіз достовірності, збіжності та ефективності розробленої методики на основі порівняння отриманих розв’язків з відомими аналітичними, чисельними результатами та експериментальними даними; Об’єкт дослідження – процеси термопружного деформування елементів конструкцій з твердих тіл. Предмет дослідження – величини параметрів напружено- деформованого стану конструкцій
в умовах термосилового навантаження. Методи досліджень – теоретичні методи: аналіз науково-технічних джерел з проблем дослідження. Емпіричні методи: аналіз проблеми чисельного моделювання фізико-механічних полів в твердих тілах; дискретизація твердих тіл чисельними методами. Експериментальні методи: тестування розробленого додатку на експериментальних задачах. Наукова новизна результатів роботи полягає у створенні на основі методу скінченних елементів нового ефективного підходу до чисельного моделювання задачі термопружності конструкцій під дією силових та температурних навантаженнях. Практичне значення отриманих результатів полягає в розробці чисельної методики і програмних засобів для визначення напружено- деформованого стану (НДС) композитних просторових конструкцій в умовах 8 термосилових навантажень з урахуванням анізотропії властивостей матеріалу. Результати магістерської роботи можуть застосовуватись при проектуванні та вдосконаленні конструкцій з композитів на відповідних промислових підприємствах машинобудівної, будівельної, транспортної, військової та інших галузей. Особистий внесок здобувача. Основні результати та положення, які становлять зміст магістерського дослідження, отримані автором самостійно. Матеріали роботи повною мірою відображені в публікаціях автора. В індивідуальній публікації, а також підготовлених у співавторстві, викладені наступні наукові результати, що належать автору: методика та алгоритм розв’язання задач деформування конструкцій з початковими напруженнями скінченного елемента в задачах термопружності [21, 6]. Апробація результатів дисертації. Матеріали дослідження доповідались на Міжнародній науковій конференції Mоdеrn wауs оf sоlvіng thе lаtеst рrоblеms іn sсіеnсе. Рrосееdіngs оf thе ХХХVІІ Іntеrnаtіоnаl Sсіеntіfіс аnd Рrасtісаl Соnfеrеnсе. (м. Vаrnа, Bulgаrіа. 2022 р.); Публікації. Основні результати дослідження викладено в 2 наукових працях, з них: 1 монографія [6], 1 стаття у виданні, що входять до закордонних міжнародних науковометричних баз [21]. Структура та обсяг роботи. Магістерська робота складається зі вступу, трьох розділів, висновків, списку використаних джерел і додатків. Загальний обсяг роботи становить 90 сторінка, у тому числі 15 рисунків, 4 таблиці, список використаних джерел із 43 найменувань на 5 сторінках, додатки на 10 сторінках. В першому розділі наведено отримані в результаті проведеного аналізу літературних джерел дані про підходи та методи розв’язання задач в умовах дії термосилових навантажень, визначено вихідні співвідношення та обрано найбільш ефективні чисельні методи їх розв’язання. 9 У другому розділі містить основні вихідні і розрахункові співвідношення методу скінченних елементів (МСЕ) для задач теплопровідності конструкцій, а також алгоритми їх розв’язання. Проведено чисельні дослідження для обґрунтування достовірності результатів. В третьому розділі наведено особливості використання МСЕ для розв’язання зв’язаних задач термопружності. Представлено структуру розробленого програмного комплексу. У додатках наведено фрагменти коду та сертифікат участі у міжнародній науково-практичній конференції. 10 РОЗДІЛ 1 АНАЛІЗ СУЧАСНОГО СТАНУ ІСНУЮЧИХ РІШЕНЬ 1.1. Проблеми та задачі термомеханіки твердих тіл В процесі експлуатації елементи конструкцій зазнають різного роду механічних і теплових навантажень, а вимоги зниження матеріаломісткості і собівартості конструкції, збільшення строків її роботи призводить до необхідності прогнозування поведінки матеріалу які неможливо описати за допомогою тільки лінійних співвідношень. В роботах В.А. Баженова [6], Я.М. Григоренка [11], О.І. Гуляра, О.С. Сахарова [13], В.В. Киричевского [19], В.Г. Пискунова [31], В.О. Толока [40], J.R. Bаrbеr [1] та ін. висвітлено термодинамічне обґрунтування основних рівнянь механіки взаємозв'язаних полів в матеріалах і елементах конструкцій, в тому числі і у неоднорідних, розвиток яких зумовлено потребам сучасної техніки внаслідок її експлуатації в умовах статичних і динамічних механічних та температурних навантажень. Як правило, елементарні термомеханічні умови об'єднують з найпростішими тепловими умовами [25]. Умови термомеханічної взаємодії з матрицею можна представити як узагальнення умов пружної взаємодії на випадок впливу температурного поля. При цьому задачі термопружності, зводяться до вирішення систем сингулярних інтегральних рівнянь щодо стрибків напружень або переміщень на лінії, яка моделює включення або тріщину, в разі, коли задані силові навантаження і температурне поле. При розв’язанні задач теплопровідності для тонкостінних систем значного спрощення можна досягти шляхом введення різних гіпотез про характер розподілу температурного поля по товщині конструкції. У роботах [34, 29] розглянуто аналіз термопружних напружень композитних матеріалів, викладено теорії і методи визначення термічних напружень. Математичне обґрунтування розрахунку деформації при впливі 11 нестаціонарного і неоднорідного по товщині пакету температурного поля багатошарових циліндричних оболонок, а також чисельного методу дискретної ортогональної прогонки запропоновано в працях [20]. Тепловий стан об'єктів, що нагріваються, можна визначати і експериментальними методами [33, 23, 16]. Проте одержана при цьому інформація є недостатньою для розрахунку полів термічного напруження, тому що має дискретний характер і температурне поле не може бути визначено в об'ємі всього тіла. До недоліків експериментальних методів слід віднести і той факт, що при вимірюванні виникає викривлення температурного поля в місцях контакту термопар з поверхнею тіла і, особливо, при спробі визначення температури у внутрішніх точках тіла. Тому у більшості випадків для аналізу теплового стану перевагу надають розрахунковим методам, порівняно з експериментальними. До класичних методів рішення задач теплопровідності належать такі методи: метод розділення змінних (метод Фур'є), метод функцій Гріна, метод теплових потенціалів і ін. Стосовно до шаруватих конструкцій при використанні методу Фур'є, спочатку визначається сукупність частинних рішень диференціального рівняння теплопровідності, а потім складається ряд з цих рішень. Сутність методу функцій джерел полягає в тому, що будь-який процес розповсюдження тепла в тілі можна представити у вигляді суми процесів вирівнювання температур, спричинених дією безлічі елементарних кількостей теплоти (джерел) розподілених в просторі і часі. Рішення для нелінійних задач в шаруватих композитах на основі вирішення у функціях Гріна розглянуто при нестаціонарних, періодичних і стаціонарних режимах теплопровідності [42]. При використанні варіаційних методів розв’язок визначаємо з умови мінімуму деякого функціоналу [7]. Розв’язок варіаційної задачі зазвичай знаходиться у вигляді, що мінімізує послідовності функцій які належать 12 області визначення еквівалентного функціоналу. Найбільш поширеними методами є методи Рітца і Канторовича [28]. Суть методу Рітца полягає в знаходженні розв’язку у вигляді ряду по базисних функціях, що мають властивість повноти і задовольняє граничним умовам задачі [28]. Невідомі коефіцієнти розкладання знаходяться шляхом мінімізації функціонала. Розв'язання задачі теплопровідності методом Канторовича, в порівнянні з методом Рітца, більш просте і більш точне. До переваг даного методу слід віднести також те, що структура наближеного розв’язку визначає апріорі тільки по одній змінній, а по іншим змінним розв’язок повністю відповідає розглянутій задачі. Застосовуються також методи Трефтца, Біо, Лейбензона, Курапка, варіаційний принцип Айонолі та ін. Докладний опис цих методів, а також використання інтегрального перетворення Лапласа і варіаційних задач наведено у працях [28, 36]. Проекційні методи або методи зважених нев'язок, стосовно до задач теплопровідності можна поділити на метод колокацій, метод найменших квадратів, інтегральний метод, метод Гальоркіна. Ці методи стосовно до задач теплопровідності докладно описано у [43, 10, 35]. Варто також зазначити і деякі інші аналітичні методи
розв’язання нелінійних задач, в яких введено поняття термічного шару: метод Швеця, метод Веймена і ін. Деякі аналітичні загальні і окремі методи розв'язання нелінійних задач теплопровідності
наведено у роботах [27]. Основною перевагою аналітичного методу є те, що розв’язок отримується у загальному вигляді. Це дозволяє визначити температурний режим в будь-який момент часу і для будь-якої точки тіла. Проте для кусково-неоднорідних (шаруватих) тіл виникають певні труднощі, пов'язані з наявністю шарів, що мають різні теплофізичні властивості. Загальний порядок вирішуваної системи диференціальних рівнянь залежить від кількості шарів. Крім того, аналітичні розв’язки можуть бути використані для 13 чітко обмеженого класу задач. Певною мірою ці задачі можна назвати класичними. У практиці проектування коло задач значно ширше. Для аналізу теплового стану реальних об'єктів, потрібне застосування чисельних методів, які дозволяють провести розрахунок об'єкту при заданих теплових діях початкових і граничних умов. Широкий розвиток комп’ютерних систем став поштовхом до розвитку ефективних чисельних методів, таких як
метод скінченних різниць (МСР), варіаційно-різницевий метод (ВРМ), метод прямих (МП), метод Монте-Карло (ММК), метод скінченних елементів (МСЕ
), метод граничних елементів (МГЕ). В першу чергу з чисельних методів необхідно відзначити МСР, який отримав значне поширення при розв’язанні задач теплопровідності. В основу МСР покладено ідею про заміну диференціальних операторів деякими різницевими аналогами. Внаслідок приходять до розв’язання систем алгебраїчних рівнянь. Питання апроксимації, збіжність, точності і стійкості МСР розглянуто в роботах [12]. Варто відзначити, що до теперішнього часу розроблено низку універсальних програм, зорієнтованих на розв’язання задач теплопровідності для однорідних тіл, в основу яких покладено МСР. Метод прямих є поєднанням МСР за однією координатою і аналітичних розв’язків за рештою. Даний метод є ефективним у окремих практичних задачах для областей простої геометричної форми. В основі методу Монте-Карло лежить моделювання статистичного експерименту на комп‘ютері і реєстрація числових характеристик, отриманих з експерименту. Серед варіаційних чисельних методів першим почали застосовувати варіаційно-різницевий метод. Головна ідея цього методу - заміна похідних під знаком інтегралу, який відповідає варіаційної постановці задачі. Цей 14 метод з успіхом застосовувався для розв'язання складних задач теплопровідності, термопружності у роботах С. Г. Міхліна [28]. Одним з найефективніших чисельних методів, що мають варіаційну основу, є МСЕ. Виникнення і розвиток МСЕ обумовлений у першу чергу роботам в галузі будівельної механіки. Застосування МСЕ дає можливість здійснювати розрахунки з урахуванням більш досконалих фізичних і математичних моделей. Для задач термопружності цей метод почав розвиватися порівняно нещодавно, на початку 60-х років минулого століття за кордоном і в 70-х роках на Україні. Великий вклад у його розвиток внесли вітчизняні та іноземні вчені В.А. Баженов [6], В.Г. Карнаухов[18], В.В. Киричевський[19], В.Г. Пискунов[32], О.С. Сахаров[19, 26], В.О. Толок[40], Л. Сегерлінд[37],та ін. Спочатку для оцінки теплового стану двовимірних тіл, почав застосовуватися трикутний скінченний елемент (СЕ) з полілінійним розподілом температури по області елементу. Надалі обґрунтованість застосування даного СЕ підтверджено в роботі [37], в тому числі для конструкцій з різко змінними геометричними і фізичними параметрами, а також в нелінійних задачах. Можливе також застосування ізопараметричних СЕ для оцінки теплового стану конструкцій [26]. Варіаційно-різницеве трактування МСЕ [36] поставило його в один ряд з іншими чисельними методами рішення рівнянь механіки та математичної фізики і дозволило довести, що рішення з використанням МСЕ збігається з порядком, визначеним порядком апроксимації інтегралу енергії. В даний час існує декілька підходів до розв’язання нестаціонарних задач теплопровідності на основі МСЕ. У разі якщо часовій аргумент є параметром варіаційних рівнянь, застосовують МСР, або метод Гальоркіна за часовою координатою [26]. Для варіаційних принципів зі згорткою проводять 15 чисельне інтегрування по часовому аргументу. У практичних застосуваннях найчастіше використовується сімейство двошарових або тришарових часових різницевих схем [15]. Для досягнення необхідної точності розв’язку задач МСЕ здійснюється послідовне згущення сітки і порівняння результатів, отриманих при відповідному збільшенні числа невідомих. При дослідженні нестаціонарних задач використовуються крокові алгоритми за параметром часу, відповідно для збільшення точності проводиться послідовне зменшенням кроків за часом. Задача теплопровідності для тіл однорідної структури на основі МСЕ розглянуто у великій кількості робіт. Нелінійні задачі розглянуто у роботі [19], а питання збіжності та стійкості розв’язків – у роботах [43]. Питання автоматизації обчислювальних процесів при розв’язанні задач МСЕ отримали висвітлення в роботах[5]. Динамічна поведінка і температура дисипативного розігріву в'язкопружного тіла обертання при осесиметричному навантаженні в нелінійній постановці розглянута в роботі [18]. Тут враховано залежність фізико-механічних властивостей матеріалу від температури, зокрема, комплексного модулю зсуву. Рішення на кожній ітерації систем лінійних рівнянь в'язкопружності і теплопровідності проводилася методом скінченних елементів. Для вирішення за МСЕ було обрано чотирикутний СЕ з квадратичною апроксимацією, що дозволило підвищити точність обчислень в області резонансу. Скінченно-елементний метод розв’язання двовимірних задач термов'язкопружності при визначенні дії температурного поля на шаруваті оболонки обертання з урахуванням нелінійних ефектів, обумовлених залежністю напружень від температури, дано в статтях [9]. На кожному ітераційному кроці вирішувалася лінеаризована задача пружності і теплопровідності. 16 Інший підхід засновано на зведенні тривимірної задачі теплопровідності до двовимірної за допомогою, будь-яких гіпотез про розподіл температури по товщині пакету шарів з подальшою реалізацією двовимірної моделі методом скінченних елементів. Таку методику було використано авторами робіт [39] при вирішенні задач теплопровідності для ортотропних та анізотропних шаруватих оболонок та пластин. Проте застосування такого підходу обмежується елементами конструкцій простої форми [31]. Існують також підходи, засновані на синтезі МСЕ і методу інтегральних перетворень [41]. Використання чисельних методів відкриває великі можливості для реалізації прикладних моделей розрахунку неоднорідних оболонок і пластин на теплові впливи. Стосовно до розрахунку шаруватих конструкцій методом скінченних елементів доцільна розробка чисельних методів з доведенням їх до програмного забезпечення, орієнтованого на широке використання в практиці проектування, хоча наявність шаруватої структури, відмінність теплофізичних характеристик вносить значні труднощі при розв’язанні задач теплопровідності. У роботі Хорошуна Л. П.[41] сформульовано модель нелінійного деформування волокнистих композитних матеріалів при фізично нелінійному деформуванні компонентів. Побудовано алгоритми для визначення ефективних деформівних властивостей волокнистих композитних матеріалів. Досліджено вплив нелінійності деформування компонентів на ефективні деформівні властивості волокнистих композитів. У роботі [1] розв’язано задачу термопружності для тонких ортотропних оболонок невід’ємної кривини при дії зосередженого джерела тепла, що рухається по поверхні оболонки. Задано лінійний розподіл температури по товщині оболонки та конвективний теплообмін за законом Ньютона з її бокових поверхонь. За допомогою інтегральних перетворень Фур’є і Лапласа отримано розв’язок в аналітичному вигляді. Досліджено вплив 17 термомеханічних властивостей матеріалу, а також параметрів теплообміну за навколишнім середовищем на напружено-деформований стан оболонки. У роботі [24] розроблено модифікований метод Максвела визначення ефективних сталих сформульовано в термінах дипольних моментів фрагмента реального композиту і еквівалентного включення. За умови врахування взаємодії між включеннями метод є строгим у сенсі, що одержана шляхом розв’язання модельної задачі оцінка ефективної сталої збігається до точного значення зі збільшенням розміру кластера. На прикладі задачі визначення теплопровідності волокнистого композиту показано, що метод забезпечує обчислення з високою точністю ефективних сталих композитів як періодичної, так і невпорядкованої мікроструктури. Задачу про визначення релаксаційних властивостей двокомпонентного матеріалу залежно від часу, концентрації армуючого компонента, типу армування, а також їх взаємодії розв'язано у роботі А.О. Камінського[17]. Тип армування визначається параметром еліпсоїда обертання, яким моделюється включення. Ефективні модулі композитного матеріалу визначено на основі релаксаційних властивостей матеріалів компонентів. Досліджено випадок, коли ізотропні матеріали обох компонентів мають в’язкопружні властивості, які можуть бути описані за допомогою двох дробово-експоненційних функцій Работнова різних порядків дробовості для характеристик об’ємного розширення й зсуву. Розв’язок в області оригіналів пропонується отримати за допомогою його дробового раціонального представлення в області перетворення. Оптимізація параметрів цього представлення й використання наведеної в роботі процедури перерахунку параметрів розв’язку в області оригіналів дозволила отримати розв’язки в компактній формі в термінах ядер релаксації. Широкий спектр методів із застосуванням імовірнісних та статистичних підходів, які дозволяють описувати механічну поведінку композитів представлено у монографії У. W. Kwоn, D. H. Аllеn, R. Tаlrеjа [2]. 18 Розглядаються різні підходи до визначення пружних сталих при наявності пошкоджень, тріщин, дифузії, температурних факторів, реологічних процесів і та ін. У тому числі наведено термопружні сталі для волокнистого композита, армованого системою n волокон. 1.2. Сучасні підходи до розв’язання задач термомеханіки у тривимірній постановці При проектуванні конструкцій, що працюють в умовах високого теплового і механічного навантаження з високими техніко-економічними показниками, виникає необхідність прогнозування фізико-механічних характеристик і показників якості композиційних матеріалів у тривимірному аспекті. Перші дослідження в області термопружності було опубліковано Пуассоном в 1828 р в роботі
«Про рівновагу та рух пружних тіл» [4]
. Трохи пізніше Дюгамель отримав рівняння термопружних коливань порожньої кулі і визначив температурні напруження всередині двошарової кулі, що знаходиться під дією змінного центрально-симетричного поля. Дослідження Дюамеля були продовжені Нейманом (1841р.) та Кельвіном (1878 г.), на основі чого було представлено рівняння лінійної термопружності, в яких рівняння стану подано лінійними залежностями між напруженнями, деформаціями і температурою. У роботах Бархардта (1873), Релея (1901), Тимошенка (1916) і ін. науковців було розроблено загальні методи розв’язання задач термопружності для двовимірного і тривимірного розподілення температур. До середини 50-их років минулого століття були вирішені деякі окремі задачі. При вирішенні просторових і плоских задач широко використовується загальні подання рівнянь теорії пружності через гармонічні, бігармонічні, квазігармонічні і інші функції. 19 У роботах Я. М. Григоренка, А. Т. Василенка, Н. Д. Панкратової [11, 12] досліджено НДС тонкостінних неоднорідних оболонок при температурних впливах, при цьому автори вводять додаткові припущення про розгляд заповнювача як тривимірного тіла, а до відносно тонких і жорстких навантажених шарів застосовуються гіпотези Кірхгофа або вони працюють як мембрани. Одним з універсальних і ефективних сучасних методів, як зазначалося раніше, при розв’язанні задач термопружності для анізотропних тіл зі складною геометрією є метод скінченних елементів. Для розв’язання просторових задач побудовано скінченні елементи різної форми - паралелепіпеди, тетраедри, трикутні призми та ін. Питанням комплексного розв’язання задач термопружності для однорідних та шаруватих тіл присвячено ряд праць. Загальні питання термопружності в механіці твердих тіл, у тому числі механіці композиційних матеріалів, знаходяться в центрі уваги дослідників. Про оптимальний опис в'язкопружних властивостей композитів, що складаються з компонентів з в'язкопружними властивостями розглянуто задачу на основі методу перетворення Лапласа – Карсона та раціональної апроксимації побудовано алгоритм зведення функцій, що описують в’язкопружні властивості матриці та наповнювача, до функцій одного класу. Отримана методика може бути використана для спрощення розв’язання задач лінійної теорії в’язкопружності. Як приклад розглянуто задачу про зведення функцій в’язкопружності до одного класу в односпрямовано армованому композиті [31]. В рамках МСЕ розв’язання задач термопружності шаруватих конструкцій можливо шляхом застосування тривимірних скінченних елементів і дискретизації всієї системи, як в плані, так і по товщині. Такий підхід застосовано у роботі [39], де розглянуто деякі задачі про термонапружений стан багатошарових пластин. 20 Розв’язання просторової задачі термопружності на основі методу Бубнова-Гальоркіна приведено в роботах [8, 13]. При цьому, досліджено термонапружений стан товстих і тришарових пластин. Надалі В. Г. Пискуновим було узагальнено метод рішення на випадок шаруватих тіл [31]. Розв’язання змішаних тривимірних задач для двошарових термопружних анізотропних плит приведено у [7]. Система рівнянь термопружності, яку записано в безмірних координатах, є сингулярно збуреною. Розв’язання здійснюється асимптотичним методом. Показано неправомірність прийняття класичної теорії пластин до розглянутого класу задач. У роботі [33] В.Н. Потураєвим проведено теоретичні та експериментальні дослідження по визначенню температурних полів і теплоутворенню у гумометалевих виробах при деформаціях стиснення, зсуву й удару в режимі циклічного навантаження. Також представлено розв’язок задачі тривимірної теплопровідності для пластини і циліндру при сталому температурному режимі, отримано вирази для інтенсивності внутрішніх джерел тепла в них при різних видах навантаження. Встановлено нелінійну залежність температури нагрівання гумових елементів від амплітуди деформації. У роботі Карнаухова В.Г. [18] розглянуто зв’язану тривимірну задачу про вимушені гармонічні коливання й дисипативний розігрів тришарової товстостінної в’язкопружної циліндричної панелі з шарнірним опиранням торців при дії на неї зовнішнього рівномірного поверхневого гармонічного тиску. Враховано взаємодію механічних і теплових полів за припущення про незалежність властивостей матеріалу від температури. Розроблено скінченно- елементний метод розв’язання динамічної задачі теплопровідності з відомим джерелом тепла. Розраховано амплітудно- й температурно-частотні характеристики, а також залежності власної частоти, коефіцієнта 21 демпфірування, максимальних амплітуди й температури розігріву від товщини демпфуючого середнього шару панелі при її коливаннях на першій моді. У роботі О.Я.Григоренка [12] на базі тривимірної теорії пружності розглянуто задачу про вільні коливання суцільного циліндру при різних граничних умовах на його торцях. Запропоновано чисельно-аналітичний підхід для розв’язання даної задачі. Рівняння теорії пружності в частинних похідних за допомогою методів сплайн-апроксимації та колокації зведено до систем звичайних диференціальних рівнянь високого порядку по радіальній координаті, які розв’язуються стійким чисельним методом дискретної ортогоналізації в поєднанні з методом покрокового пошуку. Наведено результати розрахунку для випадку трансверсально-ізотропного та неоднорідного матеріалів циліндра та деяких видів крайових умов на його торцях.У роботі [12] проведено дослідження напружено- деформованого стану пологої прямокутної в плані ортотропної оболонки зі змінною товщиною в двох координатних напрямках при різних граничних умовах на краях в уточненій постановці. Для розв’язання двовимірних крайових задач використано чисельно-аналітичний підхід, який базується на застосуванні сплайн-апроксимації та методу дискретної ортогоналізації. На теперішній час для ефективного пошуку розв’язку нестаціонарних задач потрібно володіти всім набором засобів аналізу, чітко представляти правила і границі їх застосування. До того ж, високі вимоги до розрахункових моделей, що закладені до сучасної нормативної бази припускають розгляд нестаціонарних навантажень різного рівня інтенсивності, які можуть діяти на конструкцію в один і той же момент часу, що в свою чергу потребує ретельного дослідження меж достовірного використання того чи іншого алгоритму. Важливим етапом в реалізації обчислювальних систем для розв’язання просторових задач динаміки є вибір оптимальних, з точки зору швидкості і складності процесів деформування, 22 алгоритмів інтегрування рівнянь руху у часі, особливістю яких в даній роботі є їх розвиток в рамках МСЕ. Відомо, що вибір того чи іншого алгоритму пов’язано із швидкістю протікання процесу деформування та характеру навантажень. При аналізі перехідних процесів найбільше розповсюдження отримав метод узагальнених координат, який припускає розклад розшукуваного рішення в ряд по формах власних коливань. В цілому, для розв’язання систем рівнянь динаміки рекомендовані прямі методи безпосереднього інтегрування по часовій координаті, залишаючи модальні методи для проблем, що пов’язані з використанням нижньої частини частотного спектру. Широке розповсюдження отримали неявні схеми інтегрування за часом (метод Ньюмарка, метод θ-Вільсона), які допускають у порівнянні з явними значно більший крок по параметру. Специфіка функцій у просторі оператора пружності для тіл із змінними вздовж направляючої параметрами жорсткості і мас. Явно виділена клітинна структура як матриці жорсткості, так і матриці мас, є підставою для використання алгоритмів, що комбінують прямі та ітераційні методи розв’язання. Поєднання ітераційних циклів, що з однієї сторони пов’язані з порушенням однорідності об’єкту вздовж направляючої, а з іншої є природними для самих алгоритмів (метод ітерацій у підпросторі) дозволяє значно зменшити обчислювальні задачі [6]. В роботі М. Mukhораdhуау [3] здійснено спробу викласти комплексний та єдиний підхід до аналізу ряду аспектів механіки композитів. Викладено різні теорії деформування конструкцій зі склопластиків, зокрема балок, пластин тощо. Розглянуто особливості розв’язання задач статики, динаміки, термопружності та теплопровідності, механіки руйнування. Представлено основи методу скінченних елементів для розв’язання цих задач та чисельні результати розрахунків. Проведений аналіз методів розв’язання задач теплопровідності і пружності показав, що для композитних конструкцій існують результати розв’язання в тривимірній постановці. Проте методика розрахунку 23 еластомерних композитних конструкцій потребує подальшого розвитку шляхом розробки моделей, що дозволяють описувати просторовий характер розподілу температурних полів. Найбільш перспективним є використання методу скінченних елементів, як одного з ефективних чисельних методів розв'язання задач термопружності. Зрозуміло, що практичне застосування тривимірних рівнянь теорії пружності для розрахунку реальних шаруватих конструкцій, за умов термосилових експлуатаційних навантажень з наперед заданими крайовими умовами, побудованих на основі застосування тривимірної схеми методу скінченних елементів, дозволяє ефективно моделювати поведінку шаруватих конструкцій складної форми. При цьому вдається уникнути спрощуючих гіпотез для моделювання композитних конструкцій. Враховуючи різні способи розповсюдження тепла у конструкціях, розроблено різні методи й математичні моделі розв'язання задач. Розповсюдження тепла може проходити за рахунок теплопровідності, конвекції та випромінювання. При першому способі розповсюдження тепла (теплопровідності), коли тепло передається через саму речовину, процес відбувається через нерівномірний розподіл температур у різних частинах тіла і супроводжується виникненням напряму потоку тепла і густини потоку тепла (віднесеного до одиниці площі кількості, що проходить через площадку тепла в одиницю часу). Поширення тепла внаслідок відносного руху часток нагрітого тіла завжди супроводжується виникненням в тілі температурних напружень і деформацій. Наявна залежність теплофізичних сталих матеріалу від температури і часу, залежність температури від часу, різних факторів зв'язку між полями напружень і деформацій, зв'язку теплових й електромагнітних полів і т.д. привели до розробки сучасних математичних моделей вирішення зв'язаних і незв'язаних задач термопружності. 24 При вирішенні задач термопружності для конструкцій зі анізотропних матеріалів застосовано різні теорії й підходи, які базуються на співвідношеннях зв'язаної задачі термовязкопружності, отриманих В.Г. Карнауховим [18], І.К. Сенченковим, Я.О. Жуком [14] та ін. Одним з найважливіших критеріїв дослідження пружних тіл з нерівномірним температурним полем є врахування залежності від температури фізичних і теплофізичних характеристик матеріалу: модулю пружності або модулю зсуву, коефіцієнту температурного розширення й коефіцієнту теплопровідності. В цьому випадку вирішується зв'язана задача термопружності. При використанні лінійної моделі температурний й термопружний стан визначаються рішенням системи, яка складається з рівняння теплопровідності, класичних рівнянь руху, рівнянь закону Гука і класичних рівнянь сумісності [19, 5] та ін. Очевидно, що всі термомеханічні процеси залежать від часу і для їх дослідження виникають задачі, що носять назву нестаціонарних (динамічних). Проте можна виділити ряд процесів, термомеханічний стан,в ході яких, хоча і змінюється в часі, але, починаючи з певного моменту часу, система приходить в стаціонарний (статичний) стан, який не залежить від часу. Моделі, що описують такі процеси, називають квазістатичними. У квазістатичних задачах термопружності інерційний член у рівняннях руху і похідні деформацій за часом в рівнянні теплопровідності не враховуються. Формулювання задачі для квазістаціонарної постановки можна представити у вигляді рівняння Біо і рівняння теплопровідності [30]: δFdv − δu − Qdu = 0 Р ∫∫∫ ∫∫∫ ∫∫ V V S , (1.1) 25 іj с (T T ) Tdv (T T ) dv − δ + β − δε = ∫∫∫ ∫∫∫ 0 0 іj ε V V іj [ ] T T dv w Tdv q h(T ) Tds = λ δ + δ + + − θ δ ∫∫∫ ∫∫∫ ∫∫ 0 ,і , j V V S (1.2) де F – вільна енергія; Р, Q– вектори об'ємних і поверхневих навантажень; u – вектор переміщень; – теплоємність при постійній деформації; с ε іj β – компоненти тензору ізотермічних пружних постійних, які визначають взаємний вплив температурного поля і поля деформацій; – компоненти тензору деформацій; ε іj іj – компоненти тензору теплопровідності; λ w – густина внутрішніх джерел теплоутворення; 0 q – тепловий потік; h – коефіцієнт теплообміну; θ – температура навколишнього середовища. Розв'язання задач для багатошарових конструкцій складної форми пов'язано з математичними труднощами, які пов'язані з описом геометричних параметрів, умов теплообміну з зовнішнім середовищем, умов з'єднання, урахуванням внутрішніх джерел тепла, рішенням систем інтегральних рівнянь [6]. Для практичної реалізації є найбільш прийнятними уточнені теорії, що ґрунтується на безперервно-структурних моделях, в яких кількість і порядок дозвільних рівнянь не залежить від числа шарів. Серед чисельно-аналітичних методів, широке застосування для розв’язання задач знаходять метод Бубнова-Гальоркіна, метод прямих, метод Релея-Рітца і інші. Аналітичні методи, засновані на безпосередньому інтегруванні диференціальних рівнянь, є прийнятними у випадках, коли ці рівняння і граничні умови досить прості. 26 У реальних процесах деформування дія навантаження може відбуватися тільки за якийсь проміжок часу. Миттєве прикладання навантаження - це ідеальний процес, якого в природі не існує. Якщо розглянути деформування в часі, то напруження в момент часу залежать від того, які напруження були в момент У деяких конструкціях цей фактор проявляється не суттєво. Проте існує великий клас конструкцій, в яких напруження попереднього навантаження істотно змінюють її жорсткісні властивості. Незважаючи на наявність великого досвіду розробки і експлуатації силових конструкцій, в даний час методи ідентифікації фізико- механічних параметрів розроблені недостатньо. Розглянуті вище аналітичні і чисельні результати розв’язків задач термопружності, показують, що останнім часом йде інтенсивна чисельна реалізація на основі МСЕ розроблюваних прикладних теорій. Однак кожна з методик має свою область застосування і не є універсальною. Найбільш ефективним є використання скінченно-елементної методики, побудованої на основі просторової постановки, оскільки в цьому випадку можна уникнути різних спрощуючих гіпотез і дозволити врахувати геометрію конструкції і просторовий характер термосилового навантаження. 1.3. Висновки до розділу На підставі аналізу та літературних джерел за темою магістерського дослідження можна зазначити, що незважаючи на велику кількість робіт, присвячених вивченню поведінки конструкцій з анізотропних матеріалів, залишаються ряд невирішених питань. Серед них можна виділити наступні: − у більшості випадків задачі вирішувалися за двовимірною постановкою; − в більшості робіт не враховувалася залежність фізичних і механічних властивостей матеріалу від навантаження; 27 Зроблені висновки дозволяють сформулювати основні завдання дослідження: 1. Розробка математичної моделі для чисельного моделювання конструкцій з еластомерів в задачах
визначення напружено-деформованого стану та температурних полів; 2. Створення чисельної методики моделювання термопружного деформування конструкцій; 3. Побудова алгоритмів чисельного моделювання термомеханічних полів в твердих тілах на основі розв’язання задачі термопружності; 4. Створення пакету прикладних програм для визначення напружено-деформованого стану.
28 РОЗДІЛ 2 МЕТОДИКА РОЗВ’ЯЗАННЯ ЗАДАЧ ТЕПЛОПРОВІДНОСТІ КОНСТРУКЦІЙ З ТВЕРДИХ ТІЛ 2.1. Скінченний елемент для розв’язання задач теплопровідності Основною метою результатів досліджень, викладених у цьому розділі, є побудова уточненої скінченно-елементної моделі розподілу температурного поля по товщині шаруватого анізотропного композитного тіла. Розглянемо застосування тривимірних скінченних елементів в задачах теплопровідності для шаруватих анізотропних конструкцій. Для скінченно-елементного формулювання задачі теплопровідності для просторових конструкцій з анізотропних матеріалів, що мають складну ���� 1 2 3 геометричну форму, використовується базисна ���� (���� , ���� , ���� ) і місцева і 1 2 3 ( ) ξ ξ ,ξ ,ξ ортогональні системи координат, зв'язані зі скінченнім елементом. Між базисними і місцевими координатами СЕ прийнято відображення, де прирости деформацій пов’язані співвідношеннями: ∂Х ∂Х ∂z ∂z * * m n m n ε = ε ε = ε іj mn , іj mn ∂z ∂z ∂Х ∂Х і j і j Для побудови матриці теплопровідності вихідне співвідношення опишемо у вигляді варіаційного рівняння Лагранжа при розгляді стаціонарної задачі теплопровідності: ∂T ∂δT (2.1) W w T dv q Tds h(T ) Tds 0, δ = δ + λ + δ − θ δ = о іj ∫ ∫ ∫ і j ∂ξ ∂ξ V S S 1 2 29 де - потужність джерела теплоутворення; w о θ - температура навколишнього середовища; λ - тензор теплопровідності; іj - тепловий потік; h - коефіцієнт тепловіддачі. Варіаційне рівняння (2.1) за допомогою дискретизації розрахункової схеми МСЕ зведемо до системи алгебраїчних рівнянь. Розглянемо лінійний скінченний елемент, властивості якого визначаються теплофізичними характеристиками. Розподіл температури за об’ємом СЕ задамо у вигляді лінійного закону апроксимації 8 1 2 3 (2.2) ) ( = ϕ ξ ξ ξ T , , T , ∑ ( ) ( ) l l 1 = l – температура у вузлових точках СЕ; де 1 2 3 3 3 1 1 2 1 2 - базисні функції для ( ) ( )( )( ) 1 1 1 + ξ ξ + ξ ξ + ξ ξ ϕ ξ ,ξ ,ξ = ( ) ( ) ( ) l l l ( l ) 4 апроксимації температури в межах елемента. Вважаємо, що між шарами матеріалу виконуються умови ідеального теплового контакту. Отже, на кордоні k-го і (k + 1)-го шару виконується умова: ( ) ( ) + + − k k 1 , = T T ( ) ( ) ( ) ( ) 1 1 k k + + + − k k (2.3) λ = λ T T . , 33 3 33 3 , Для скінченних елементів, що моделюють процес поширення тепла в двошаровому матеріалі, поле температур апроксимується за допомогою відповідних базисних лінійних функцій, а в напрямку, перпендикулярному 30 границі шарів - кусково-лінійною, що забезпечує тотожність теплових потоків в напрямку спільної границі шарів: −1 3 ξ + 1 h h і і 1 2 , φ = 2 −і + −
37%
33%
20%
12%
7%
7%
6%
4%
3%
3%
2%
2%
2%
2%
2%
1%
1%
1%
1%
0,9%
0,9%
0,9%
0,8%
0,7%
0,7%
0,1%
0%
1 + ( ) 3 1 1 2 λ λ λ 33 33 33 −1 3 h ξ − h + 1 h h і k 1 1 1 2 , і, k=1,2. 2 1 φ = −k + − + + ( ) 3 1 2 1 2 λ λ λ λ 33 33 33 33 Для 1-го шару −1 3 h h 1 1 ξ + 1 2 3 1 2 2 1 ( , , ) ( ϕ ξ ξ 1 )( ξ = 1 ) 1 − ξ − ξ − + , (1) 2 1 1 4 λ λ λ 33 33 33 −1 3 h h 1 1 ξ + 1 2 3 1 2 2 1 ( , , ) (1 )(1 ) 1 + ϕ ξ ξ ξ = + ξ − ξ − , (2) 2 1 1 4 λ λ λ 33 33 33 −1 3 h h 1 1 ξ + 1 2 3 1 2 1 2 ( , , ) (1 )(1 ) 1 + ϕ ξ ξ ξ = − ξ + ξ − , (3) 1 1 2 4 λ λ λ 33 33 33 −1 3 h h ξ + 1 1 3 1 2 1 2 2 1 + + ξ + ξ − ϕ ξ ξ ξ = 1 )(1 ) 1 ( , , ) ( , (4) 1 2 1 4 λ λ λ 33 33 33 1 − 3 h h 1 1 + ξ 1 2 3 1 2 2 1 , 1 )(1 ) ( , , ) ( + ξ ϕ ξ ξ ξ = − − ξ (5) 2 1 1 4 λ λ λ 33 33 33 −1 3 h h ξ + 1 1 1 2 3 1 2 1 2 , ϕ ξ ξ ξ = + ξ − ξ + ( , , ) (1 )(1 ) (6) 1 1 2 4 λ λ λ 33 33 33 −1 3 h h ξ + 1 1 1 2 3 1 2 1 2 , ϕ ξ ξ ξ = − ξ + ξ + ( , , ) (1 )(1 ) (7) 1 1 2 4 λ λ λ 33 33 33 1 − 3 h h 1 1 ξ + 1 2 3 1 2 1 2 . (2.4) ( , , ) (1 )(1 ) + ϕ ξ ξ ξ = + ξ + ξ (8) 1 1 2 4 λ λ λ 33 33 33 31 Для 2-го шару −1 3 h h 1 h h ξ − + 1 1 2 3 1 2 1 1 1 2 ( , , ) (1 )(1 ) 1 + ϕ ξ ξ ξ = − ξ − ξ − + , (1) 1 2 1 2 4 λ λ λ λ 33 33 33 33 −1 3 h h 1 h h ξ − + 1 1 2 3 1 2 1 1 1 2 ( , , ) (1 )(1 ) 1 + ϕ ξ ξ ξ = + ξ − ξ − + , (2) 1 2 1 2 4 λ λ λ λ 33 33 33 33 −1 3 h h 1 h h ξ − + 1 1 2 3 1 2 2 1 1 1 ( , , ) (1 )(1 ) 1 + ϕ ξ ξ ξ = − ξ + ξ − + , (3) 2 1 1 2 4 λ λ λ λ 33 33 33 33 −1 3 − + ξ 1 h h h h 1 1 2 2 3 1 2 1 1 1 + + − ξ + ξ = + ξ ξ ϕ ξ ) 1 1 )(1 ( ) , , ( , (4) 2 1 2 1 4 λ λ λ λ 33 33 33 33 −1 3 1 h ξ − h + h h 1 1 2 3 1 2 1 1 1 2 , ( , , ) (1 )(1 ) + ϕ ξ ξ ξ = − ξ − ξ + (5) 1 2 1 2 4 λ λ λ λ 33 33 33 33 −1 3 h h h h ξ − + 1 1 1 2 3 1 2 1 1 1 2 , ϕ ξ ξ ξ = + ξ − ξ + + ( , , ) (1 )(1 ) (6) 1 2 1 2 4 λ λ λ λ 33 33 33 33 −1 3 1 h ξ − h + h h 1 1 2 3 1 2 1 1 1 2 , ( , , ) (1 )(1 ) + ϕ ξ ξ ξ = − ξ + ξ + (7) 1 2 1 2 4 λ λ λ λ 33 33 33 33 −1 3 1 h ξ − h + h h 1 1 2 3 1 2 1 1 1 2 . (2.5) ( , , ) (1 )(1 ) + ϕ ξ ξ ξ = + ξ + ξ + (8) 1 2 1 2 4 λ λ λ λ 33 33 33 33 Для скінченних елементів, що моделюють процес поширення тепла в тришаровому матеріалі, поле температур апроксимується за допомогою базисних функцій вигляду: 32 Для 1-го шару: −1 3 h h h 1 1 ξ + 1 2 3 1 2 3 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = − ξ − ξ − , (1) 1 1 2 3 4 λ λ λ λ 33 33 33 33 −1 3 h h h 1 1 ξ + 1 2 3 1 2 3 2 1 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = + ξ − ξ − , (2) 2 3 1 1 4 λ λ λ λ 33 33 33 33 −1 3 h h h 1 ξ + 1 1 2 3 1 2 3 1 2 + + ξ , ξ , ξ ) = (1 − ξ )(1 + ξ ) 1 − ϕ ( , (3) 1 1 2 3 4 λ λ λ λ 33 33 33 33 −1 3 h h h 1 1 ξ + 1 2 3 1 2 3 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = + ξ + ξ − , (4) 1 1 2 3 4 λ λ λ λ 33 33 33 33 −1 3 h h h 1 1 ξ + 1 2 3 1 2 3 1 2 , ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = − ξ − ξ (5) 1 1 2 3 4
λ λ λ λ 33 33 33 33 −1 3 h h h 1 1 ξ +
37%
33%
20%
12%
7%
7%
6%
4%
3%
3%
2%
2%
2%
2%
2%
1%
1%
1%
1%
0,9%
0,9%
0,9%
0,8%
0,7%
0,7%
2 3 1 2 1 3 1 2 , + + ( , , ) (1 )(1 ) ϕ ξ ξ ξ = + ξ − ξ (6) 1 1 2 3 4 λ λ λ λ 33 33 33 33 −1 3 h h h 1 1 ξ + 1 2 3 1 2 3 1 2 , + + ( , , ) (1 )(1 ) ϕ ξ ξ ξ = − ξ + ξ (7) 1 1 2 3 4 λ λ λ λ 33 33 33 33 −1 3 h h h 1 1 ξ + 1 2 3 1 2 3 1 2 . (2.6) + + ( , , ) (1 )(1 ) ϕ ξ ξ ξ = + ξ + ξ (8) 1 1 2 3 4 λ λ λ λ 33 33 33 33 Для 2-го шару: −1 3 h h h 1 h h ξ − + 1 1 2 3 1 2 3 1 1 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = − ξ − ξ − + , (1) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 33 −1 3 h h h 1 h h ξ − + 1 1 2 3 1 2 3 2 1 1 1 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = + ξ − ξ − + , (2) 2 3 1 1 2 4 λ λ λ λ λ 33 33 33 33 33 −1 3 h h h 1 h h ξ − + 1 1 2 3 1 2 3 1 1 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = − ξ + ξ − + , (3) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 −1 3 h h h 1 h h ξ − + 1 1 2 3 1 2 3 1 1 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = + ξ + ξ − + , (4) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 −1 3 h 1 h ξ − h + h h 1 1 2 3 1 2 3 1 1 1 2 , (1 )(1 ) ( , , ) + + ϕ ξ ξ ξ = − ξ − ξ + (5) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 −1 3 h 1 h ξ − h + h h 1 1 2 3 1 2 3 1 1 1 2 , ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = + ξ − ξ + (6) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 −1 3 h 1 h ξ − h + h h 1 1 2 3 1 2 3 1 1 1 2 , ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = − ξ + ξ + (7) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 1 − 3 h 1 h ξ − h + h h 1 1 2 3 1 2 3 1 1 1 2 . (2.7) ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = + ξ + ξ + (8) 1 2 1 2 3 4 λ λ λ λ λ 33 33 33 33 33 Для 3-го шару: −1 3 1 h h ξ − + h h h h 1 1 2 3 1 2 3 3 1 2 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = − ξ − ξ − + + , (1) 2 3 1 3 1 2 4 λ λ λ λ λ λ 33 33 33 33 33 33 −1 3 1 h h ξ − + h h h h 1 1 2 3 1 2 3 3 1 2 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = + ξ − ξ − + + , (2) 2 3 1 3 1 2 4 λ λ λ λ λ λ
33 33 33 33 33 33 −1 3 1 h h ξ − + h h h h
37%
33%
20%
12%
7%
7%
6%
4%
3%
3%
2%
2%
2%
1%
0,9%
0,9%
0,9%
0,5%
1 1 2 3 1 2 3 3 1 2 1 2 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = − ξ + ξ − + + , (3) 2 3 1 2 3 1 4 λ λ λ λ λ λ 33 33 33 33 33 33 34 −1 3 1 h h ξ − + h h h h 1 1 2 3 1 2 3 3 2 1 2 1 ( , , ) (1 )(1 ) 1 + + ϕ ξ ξ ξ = + ξ + ξ − + + , (4) 2 3 1 2 3 1 4 λ λ λ λ λ λ 33 33 33 33 33 33 −1 3 1 ξ − + h h h h h h 1 1 2 3 1 2 3 3 2 1 2 1 , ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = − ξ − ξ + + (5) 1 2 3 1 2 3 4 λ λ λ λ λ λ 33 33 33 33 33 33 1 − 3 1 ξ − + h h h h h h 1 1 2 3 1 2 3 3 1 2 1 2 , ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = + ξ − ξ + + (6) 1 2 3 1 2 3 4 λ λ λ λ λ λ 33 33 33 33 33 33 1 − 3 1 ξ − + h h h h h h 1 1 2 3 1 2 3 3 1 2 1 2 , ( , , ) (1 )(1 ) + + ϕ ξ ξ ξ = − ξ + ξ + + (7) 1 2 3 1 2 3 4 λ λ λ λ λ λ 33 33 33 33 33 33 1 − 3 h h ξ − + 1 h h h h 1 1 2 3 1 2 3 3 1 2 1 2 ( ) ( )( ) , , . ϕ ξ ξ ξ = + ξ + ξ + + + + 1 1 (2.8) ( ) 8 3 3 1 2 1 2 4 λ λ λ λ λ λ
33 33 33 33 33 33 Вираз (2.2) запишемо в матричному вигляді: T { } { } = φ T T . (2.9) Апроксимуємо функцію температури у вигляді розкладання по степеневим функціям: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 000 100 100 010 010 001 001 110 110 = χ + χ ψ + χ ψ + χ ψ + χ ψ T ( ) ( ) ( ) ( ) ( ) ( ) 101 101 011 011 111 111 + χ ψ + χ ψ + χ ψ ; (2.10) або 111 рqr ) ( рqr ) ( T = Х ψ (2.11) ∑ рqr Функції форми можна представити через степеневі функції: 35 [ ]{ } { } φ = B ψ , (2.12) - матриця перетворення. де Похідні температури за координатами обчислюються за формулами: T δ ( ) ( ) ( ) ( ) ( ) ( ) ( ) 100 110 010 101 001 111 011 ; = χ + χ ψ + χ ψ + χ ψ 1 ∂ξ T δ ) ( ) ( ) ( ) ( ( ) ( ) ) (010 110 100 011 001 111 101 ; = χ + χ ψ + χ ψ + χ ψ 2 ∂ξ δT ( ) ) ) ( ( ( ) ) ) ( ) ( (001 101 100 011 010 111 110 (2.13) = χ + χ ψ + χ ψ + χ ψ . 3 ∂ξ Представимо вирази (2.12) і (2.13) в матричному вигляді і запишемо їх варіації: T { } { } δ = T δχ ψ ; T } { { } δ = δχ ψ T . (2.14) ,і ,і Згідно (2.12) маємо: T T [ ] [ ] { } { } { } { } = ψ = ψ T T B ;T T B ; ,і ,і T T T T [ ] [ ] { } { } { } { } T B T ;T B T ; δ = ψ δ = ψ δ (2.15) ,і ,і Припустимо розподіл температури на поверхні скінченного елемента, що знаходиться на кордоні області, функцією: 36 1 1 ( ) ( ) р q р q ( ) st s t t s (2.16) = φ . T T ∑ ∑ =0 =0 р q s t З урахуванням температура в матричній формі на поверхні буде дорівнювати: T T ( ) ( st ) ( st ) ( st ) st , (2.17) { } { } [ ] = ψ T B T , T T ( ) ( st ) ( st ) ( st ) st . (2.18) { } { } [ ] δT = ψ B δ T , Підставимо вираз (2.17), (2.18) у рівняння (2.1): 1 T T T T 1 2 3 δ = δ ψ λ ψ ξ ξ ξ + W T B B T gd d d } { { } ] [ [ ] } { { } і іj
37%
33%
20%
12%
7%
7%
6%
0,3%
0,1%
j , , ∫∫∫ − 1 1 T 1 2 3 ] [ { } { } w T B g d d d + δ ψ ξ ξ ξ + ∫ ∫ ∫ 0 1 − 1 s t T nn 1 1 [ ]{ } { } + δ ψ ξ ξ + q T B g Є d d ∫∫ s t s t ns t ` 1 1 1 1 1 1 1 − 1 s t T T T nn 2 2 [ ] [ ] { }{ } { } { } + δ ψ ψ ξ ξ − h T B B T g Є d d ∫∫ s t s t s t s t ns t ` ` 2 2 2 2 2 2 2 2 2 2 1 − 1 s t T nn 2 2 [ ]{ } { } d , h T B g Є d − θδ ψ ξ ξ =0 ∫∫ ns t s t s t ` 2 2 2 2 2 2 1 − ( ) s t 2 2
[ ] [ ] { } { }{ } { }{ } { } { } { } δ + δ + δ + δ = 0 T H T T H T T Р T Р , (2.19) де g– метричний тензор; Є – символ Веблена. nts Матриця теплопровідності [H] має вигляд: 37 1 T T 1 2 3 { } { } ] [ ] [ [ ] H = B ψ λ ψ B g dξ dξ dξ (2.20) ∫∫ ∫ ,і іj ,і −1 ���� ���� 2 2 ] [ ] [ є доповненням до матриці ���� , обумовленої Матриця ���� граничними умовами 3 роду на поверхні СЕ. { } Вектор еквівалентного навантаження ���� , обумовлений внутрішнім джерелом теплоутворення, надається виразом: 1 1 2 3 [ ]{ } { }= ψ ξ ξ ξ Р w B g d d d . (2.21) ∫∫∫ 0 1 − , обумовлений тепловими Вектор еквівалентного навантаження { } S потоками і температурою на поверхні СЕ, запишемо у матричному вигляді: 1 s t nn 1 1 ] [ } { } {ψ ξ ξ S q B g Є d d = −
s t s t ns t ∫ ∫ 1 1 1 1 1 1 1 (2.22) 1 s t nn 2 2 ] [ } { θ ψ ξ ξ h B g Є d d . − s t s t ns t ∫ ∫ 2 2 2 2 2 2 1 Якщо вважати, що варіація функції поля (2.19) дорівнює нулю, то і коефіцієнти при варіаціях температур також дорівнюватимуть нулю: ( s t ) 2 2
[ ] [ ] { } { } { } { }+ + + = 0 H T H T Р S . (2.23) Тоді система розв`язувальних рівнянь стаціонарної задачі теплопровідності має вигляд: 38 [ ]{ } { } ���� ���� = − ���� , (2.24) { } { } { } де ���� = ���� + ���� - вектор еквівалентного навантаження; { } ���� - поле температур. При обчисленні матриці теплопровідності необхідно виконати перетворення тензора теплопровідності, заданого в системі координат за армування анізотропного шару, в місцеву систему координат ���� ���� формулою [36]: * m n , (2.25) λ = λ а а іj mn і j ���� де ���� - тензор повороту системи координат, який пов’язаний з тензором ���� і d dх і і m і наступним співвідношенням: перетворення координат , а = d = m m dξ g m mn g ( - компоненти метричного тензора). mn При обчисленні матриці теплопровідності для скінченних елементів з кусково-лінійним законом апроксимації переміщень по об’єму елемента для забезпечення достатньої точності необхідно процедуру інтегрування проводити для кожного шару окремо. У такому випадку процедура формування системи розв’язувальних рівнянь для пакету шарів конструкції в цілому може бути зведена до суперелементного підходу, що є виправданим при розрахунку багатошарових конструкцій. При побудові системи рівнянь для суперелементної схеми необхідно перейти до глобальної нумерації вузлів кожного шару. Після компонування отримуємо систему, що складається з 12(n+1) рівнянь (n – кількість шарів). Для визначення переміщень пакету необхідно спочатку визначити переміщення на зовнішній поверхні, а саме вузлів
20%
7%
3%
1%
0,9%
0,9%
0,9%
1, 2, 3, 4, (4n+1), (4n+2), (4n+3), (4n+4). 39 Поле вузлових температур можна виразити у вигляді одновимірного масиву: T T T = Q Q Q , { } { } 1 2 3 T T Q = T T T T , { } { } 1 1 2 3 4 T T Q = T T T T
, } { { } 2 5 6 4n −1 4n T T Q =
T T T T . { } { } 3 4n +1 4n + 2 4n +3 4n + 4 Вектор еквівалентного навантаження також можна представити у вигляді сукупності векторів: T T , R = Р Р Р { } { } 1 2 3 T T Р = R R R R , { } { } 1 1 2 3 4 T T Р R R R R = … , { } { } 2 5 6 4n 1 4n − T T Р = R R R R . (2.
26) { } { } 1 4n +1 4n + 2 4n +3 4n + 4 Тоді головна система рівнянь записується наступним чином: H H H Q Р 11 12 13 1 1 H H H Q Р = . (2.27) 21 22 23 2 2 H H H Q Р 31 32 33 3 3 Поле температур вузлів обчислюється через температуру Q 1 −1 −1 Q = H Р H H− H [ ] 2 22 2 22 21 23 Q 3 40 Q 1 -1 -1 H Q + H H Р - H H H + H Q = Р ; [ ] 11 1 12 22 2 22 21 23 13 3 1 Q 3 Q 1 -1 -1 H Q + H H Р - H H H + H Q = Р ; (2.28) [ ] 31 1 32 22 2 22 21 23 33 3 3 Q 3 Таким чином отримуємо систему рівнянь для скінченного суперелемента: = H Q F , (2.29) { } } { H G H H 12 11 13 H = - - матриця теплопровідності суперелемента; де H H H G 31 33 32 -1 ; G = H H H [ ] 22 21 23 Q 1 = Q – вектор вузлових температур на зовнішніх гранях } { Q 3 16 суперелемента; Р K S 1 12 F = - { } – вектор еквівалентного навантаження на зовнішніх Р K S 3 32 -1 гранях суперелемента. . S = H Р 22 2 Процедуру обернення матриці, враховуючи її симетричність, можна виконати за методом Халецького. 2.2. Нестаціонарні задачі теплопровідності Зв'язана задача термопружності в загальному випадку формулюється для динамічного навантаження. В цьому випадку процес поширення тепла в залежності від часу, тобто слід вирішувати задачу нестаціонарної теплопровідності. Застосування аналітичних методів для таких завдань пов'язано з певними труднощами, зумовлені різними причинами, наприклад, 41 досить складною геометричною формою конструктивного елементу або складним характером теплового навантаження. Розглянемо алгоритм чисельного розв’язання нестаціонарної задачі теплопровідності. Розподіл температур по товщині пакета анізотропних шарів в цьому випадку можна вважати лінійним на кожному шарі. В цьому випадку описаний вище СЕ для розв’язання стаціонарного рівняння теплопровідності можна використовувати і для нестаціонарної теплопровідності. Рівняння нестаціонарної теплопровідності для анізотропного тіла можна представити у вигляді [34]: t t t t
2 2 2 2 іj ρ = λ + + − − θ сT Tdvdt T T dvdt w Tdvdt q h T dsdt , (2.30) ) ( ,t ,і , j 0 ∫∫∫∫ ∫∫∫∫ ∫∫∫∫ ∫∫∫ t V t V t V t S 1 1 1 1
де ρ - питома вага матеріалу; с - теплоємність; іj λ - тензор теплопровідності; - потужність внутрішніх джерел тепла; w 0 q - інтенсивність теплових потоків; h - коефіцієнт теплопередачі; θ - температура навколишнього середовища. До цього рівняння необхідно додати початкові та граничні умови. Розв’язання рішення задачі нестаціонарної теплопровідності вимагає завдання температури у всьому тілі на початковий момент часу t = 0 і ( ) = T T t ,z , 0 42 а також граничних умов. Крайові умови для задачі теплопровідності можуть бути: першого роду (задано розподіл температури по поверхні S тіла у довільний момент часу): і ( ) T = ϕ t ,z S другого роду (задана густина теплового потоку в кожній точці поверхні тіла) k ( ) T t , z ∂ k ( ) ( ) q t , z T ; = −λ q ∂n третього роду (задано конвективний теплообмін між поверхнею тіла і навколишнім середовищем): k ( ) ∂T t ,z k k [ ( ) ( )] ( ) ( ) − θ = −λ h t T t ,z t ,z T . T q ∂n ( ) h t де – температура зовнішнього середовища; – коефіцієнт теплообміну. θ Т Застосовуючи скінченно-елементну дискретизацію за координатами і часу рівняння (2.30) можна привести до системи алгебраїчних рівнянь. Функцію температури можна представити у вигляді: T ( х, t ) = f ( х) g (t ) (2.31) Для скінченного елемента, що моделює поведінку шаруватого елемента конструкції, функція координат може бути представлена у вигляді (2.2) Апроксимацію температури всередині СЕ в залежності від часу приймаємо: за лінійним законом: 43 g t = γ + γ
37%
33%
20%
12%
7%
6%
4%
3%
3%
2%
0,5%
0,4%
t ; ( ) 1 1 2 за квадратичним: 2 g t = γ + γ t + γ t ; ) ( 2 1 2 3 за кубічним законом: 2 3 = γ + γ + γ + γ g t t t t . ( ) 3 1 2 3 4 Тут – коефіцієнти розкладання, що визначаються через задані часові α і інтервали. Прийнята апроксимація дозволяє при обчисленні коефіцієнтів матриці теплопровідності інтеграли, що залежать від часу, обчислювати окремо: t 2 А = g t dt ( ) і і ∫ t 1 γ 2 2 2 ; А t t t t γ = + − ( ) ( ) 1 1 2 1 2 1 − 2 γ γ 2 2 3 3 3 2 ; А t t t t t t γ = + − + − ( ) ( ) ( ) 2 1 2 1 2 1 2 1 − 2 3 γ γ γ 2 2 3 3 4 4 3 2 4 . γ А = +t t − t + t − t + t − t t ( ) ( ) ( ) ( ) 3 1 2− 1 2 1 2 1 2 1 2 3 4
B Коефіцієнти можна обчислити з доданка, що входить в формулу: і ) ( t 2 а b ( ) ( ) B = T T dt , ,t і ∫ ( ) t 1 ( ) а ( ) b а де ; - похідна від температури за часом в стані і значення T T ,t . температури в стані b Для лінійної апроксимації за часом маємо γ β
37%
33%
20%
12%
7%
6%
2%
2%
2%
1%
0,5%
0,4%
0,1%
2 2 2 2 B =γ β t t + t − t ( ) ) ( 1 2 1 2− 1 2 1 2 44 Для квадратичної або кубічної функції g коефіцієнти B обчислюються і і аналогічно 1 1 2 2 3 3 =γ β + γ β + γ β − + γ β + γ β − + 2 2 B t t t t t t ( ) ( ) ( ) ( ) ( ) 2 2 1 2 1 3 1 2 2 2 1 2 3 3 2 2 1 − 2 3 1 4 4 + γ β − t t ; ( ) 3 3 2 1 2 1 1 2 2 3 3 B t t 2 t t 2 3 t t =γ β + γ β + γ β − + γ β + γ β + γ β − + ( ) ( ) ( ) ( ) ( ) 3 2 1 2− 1 3 1 2 2 2 1 2 3 3 2 4 4 2 1 2 3 1 1 4 4 5 5 + γ β + 3γ β + 2γ β t − t + 2γ β + 3γ β t − t + ( ) ( ) ( ) ( ) 2 4 4 2 3 3 2 1
3 4 4 3 2 1 4 5 1 6 6 + γ β t − t . (2.32) ( ) 4 4 2 1 6 - коефіцієнти розкладання, що визначаються через прийняті проміжки тут β і . часу, аналогічно γ і D Коефіцієнти знаходимо з виразу і ( ) t 2 а в ) ( ) ( D = T T dt ; і ( ) ∫
37%
33%
20%
12%
7%
7%
6%
3%
3%
2%
2%
2%
1%
1%
0,8%
0,7%
0,5%
0,4%
0,3%
0,1%
t 1 1 1 2 2 3 3 ; D = γ β t t + γ β + γ β t − t + γ β t − t ( ) ( ) ( ) ( ) − 1 1 1 2 1 1 2 2 1 2 1 2 2 2 1 2 3 1 1 2 2 3 3 D t t t t t t = γ β + γ β + γ β − + γ β + γ β + γ β − + ( ) ( ) ( ) ( ) ( ) 2 1 1 2 1 1 2 2 1 2 1 1 3 3 1 2 2 2 1 − 2 3 1 1 4 4 5 5 ; t t t t + γ β + γ β − + γ β − ( ) ( ) ( ) 2 3 3 2 2 1 3 3 2 1 4 5 1 1 2 2 3 3
D = γ β t t + γ β + γ β t − t + γ β + γ β + γ β t − t + ( ) ( ) ( ) ( ) ( )
37%
33%
20%
7%
3%
1%
0,1%
3 1 1 2− 1 1 2 2 1 2 1 1 3 3 1 2 2 2 1 2 3 1 4 4 + γ β + γ β + γ β + γ β t − t + ( ) ( ) 1 4 4 1 2 3 3 2 2 1 4
1 1 5 5 6 6 + γ β + γ β + γ β t − t + γ β + γ β t − t + ( ) ( ) ( ) ( ) 2 4 4 2 3 3 2 1 3 4 4 3 2 1 5 6 1 7 7 + γ β t − t . (2.33) ( ) 4 4 2 1 7 45 Матрицю теплопровідності отримуємо з формули (2.30) з урахуванням (2.32) для відповідного закону апроксимації: 1 і а b ( ) ( ) ( ) kl kl 1 2 3 = λ H D g T T gdх dх dх , і ,k ,l ∫ ∫ ∫ ( ) аb ( ) −1 а b ) ) ( ( T T де , - похідні від температури по місцевим координатам СЕ в ,k ,l станах а та b. Аналогічно знаходимо матриці теплоємності 1 і а b ( ) ( ) ( ) 1 2 3 G = B сT T ρ gdх dх dх , t , і ∫ ∫ ∫ ( ) аb ( ) −1 граничних умов 1 і а b ( ) ( ) ( ) 1 2 S = D hT T g dх dх , s s s і ∫ ∫ ( ) аb ( ) −1 а також відповідний вектор правої частини системи рівнянь: 1 1 і ( ) 1 2 3 1 2 . А А T gdх dх dх q h T g dх dх = ω + − θ ( ) 0 s s s і а а ∫ ∫ ∫ ∫ ∫ ( ) ( ) ( ) аа ( ) − − 1 1 тут - метричний тензор поверхні; - інтенсивність теплових потоків, g q s нормальних до поверхні тіла. Розв’язувальна система рівняння нестаціонарної теплопровідності буде мати вигляд: і і і і ( ) ( ) ( ) ( ) H + G + S + А = 0 . (2.34) аb аb аb аа ( ) ( ) ( ) ( ) Нехтуємо в ній залежністю фізико-механічних властивостей еластомеру від часу і температури. Розмір загальної матриці 46 теплопровідності залежить від законів апроксимації функції температур і координати за часом. Коефіцієнти матриць системи рівнянь і векторів правих частин визначається чисельно на основі квадратурної формули Гауса за поліномами D і А і B і Лежандра. Попереднє обчислення величин , , , звільняє від ( ) ( ) ( ) необхідності виконувати чисельне інтегрування по координаті . t Застосування просторово-часового апроксимації з використанням неявної різницевої схеми Кранка – Николсона приводить, з одного боку, до зростання порядку системи рівнянь, а з іншого – збільшує точність наближених обчислень при значно більших кроках за часом у порівнянні з явною різницевою схемою. Розв’язання задачі нестаціонарної теплопровідності зводиться до обчислювальної процедури, яка передбачає покрокове визначення температури при відповідному часовому інтервалі: 1 1 1 [ ] [ ] [ ] [ ] { } { } + H G T = G T − R , (2.35) n n n +1 ∆t ∆t ∆t де H, G – матриці теплопровідності і теплоємності; ∆t - крок за часом; Процес розв’язання системи рівнянь триває до тих пір, поки виконається умова у вигляді: T − T ≈ ε (2.36) n+1 n де ε – точність розрахунків. 2.3. Дослідження збіжності рішення Для обґрунтування достовірності пропонованого підходу розглянемо розв’язання задач стаціонарної теплопровідності для шаруватих конструкцій. 47 Задача 2.1. Розглянемо розв’язання задачі стаціонарної теплопровідності тришарової панелі. Панель прямокутна: розміри панеліа =3,5м,а =2,9м, а =0,35м. Матеріал шарів ізотропний. Товщина шарів 1 2 3 k ( ) 0,17;0,1;0,08м. Конструкція являє собою стінову панель житлового h = будинку. Температура середовища на зовнішній поверхні θ =318K, на 1 =293K. Коефіцієнти теплообміну на зовнішній і внутрішній поверхні θ 2 2 2 8,7Вт/(м ∙K), 23Вт/(м ∙K) відповідно. внутрішній поверхнях h = h = 1 2 k k k ( ) ( ) ( ) λ λ= λ= = 0,52; 0,12; 0,52. Бічні Коефіцієнти теплопровідності 11 22 33 поверхні панелі теплоізольовані. Результати розрахунку при різних сітках скінченних елементів представлено в таблиці 2.3 і на рис.
2.1, 2.2. Таблиця 2.1.- Температура на поверхнях тришарової панелі 0,1а ; 0,2а ; 0,3а ; 0,4а ; 0,5а ; 1 1 1 1 1 Координати точок поверхні 0,1а 0,2а 0,3а 0,4а 0,5а 2 2 2 2 2
Температура на поверхні 293,82 295,43 295,48 295,49 29,49 панелі МССЕ (4х11х11) 313,41 316,87 317,05 317,06 317,06 294,13 295,41 295,48 295,49 296,49 МССЕ (2х16х16) 316,73 317,01 317,06 317,06 317,06 294,17 295,41 295,49 295,49 296,49 МССЕ (2х21х21) 316,74 317,04 317,06 317,06 317,06 294,4 294,9 294,9 294,9 294,9 Розв’язок [39] (сітка 40×40) 317,0 317,2 317,3 317,3 317,3 291,69 294,29 294,65 294,70 294,71 Розв’язок [8] 316,06 317,28 317,35 317,35 317,35 48 Рис. 2.2. Розподіл температур Рис. 2.1. Розподіл температур ΔТ при сітці 4х21х21 ΔТ при сітці 4х16х16 Результати задовільно збігаються з аналітичним рішеннями, наведеними в роботах [8, 39], отриманими на основі двовимірної моделі. Перевірка збіжності результатів проведена при різних сітках скінченних елементів. При збільшенні числа скінченних елементів спостерігається стійка збіжність результатів, що свідчить про достовірність отриманих рішень. Задача 2.2. Розглянемо п’ятишарову квадратну плиту 2 × 2 м, на бічних поверхнях якої підтримується нульова температура. На нижній поверхні плити має місце конвективний теплообмін з навколишнім середовищем, 2 коефіцієнт теплообміну 10Вт/(м ∙K), температура навколишнього ◦ 0 15 С . До верхньої поверхні підводиться теплової потік середовища 2 інтенсивністю 200Вт/м . Кожен шар пластини - односпрямований епоксидний вуглепластик. Коефіцієнти теплопровідності відповідно рівні 14,6Вт/(м∙K), 0,93Вт/(м∙K). Шари орієнтовані по черзі під кутом 0° та 90°. Товщина шарів 0,17а; 0,25а; 0,16а; 0,25а; 0,17а м, де 0,1м. Результати розв’язання даної задачі, отримані при різних сітках СЕ, наведено в таблиці 2.2 і на рис. 2.3. 2.3. 49 Таблиця 2.2 - Температура в центральних точках на поверхнях п'ятишарової плити Сітка СЕ Точний Температура на розв’язок [38] поверхні панелі 6×4×4 6×6×6 6×10×10 31,314 30,720 30,043 29,906 50,540 49,749 48,692 48,457 Рис. 2.3. Розподіл температур Рис. 2.4. Розподіл температур на поверхні плити по товщині плити Задача 2.3. П'ятишарова трапецевидна панель під дією теплового навантаження. Розміри плити в плані: а=150см, b= 81см, с=15см. До верхньої 2 поверхні =12мм підводиться теплової потік =2,24кВт/м . На нижній поверхні ( = 0) і бічних гранях панелі –2-3, 3-4, 1-4 має місце конвективний =153 К), грань 1- 2 теплообмін з навколишнім середовищем ( теплоізольована. Шари 2, 4 - трансверсально-ізотропні, а 1, 3, 5 - ортотропні. Головні осі ортотропії шарів 1, 3, 5 розташовані під кутом 30° до прийнятої системі координат , отже, в цій системі характеристики матеріалу х ' Ох ' 1 2 шару відповідають випадку прямолінійної анізотропії. Теплофізичні характеристики пакета шарів (k=1…5) наступні: =0,12; 16,24; 0,12; 16,24; 0,12 Вт/(м∙К); =0,54; 16,24; 0,54; 16,24; 0,54 Вт/(м∙К); 50 =0,048; 0,83; 0,048; 0,83; 0,048 Вт/(м∙К); = =0,1785; 0; 0,1785; 0; 0,1785 Вт/(м∙К). Коефіцієнти теплообміну по бічних гранях =1500; 1700; 1500; 1700; 1500Вт/(м∙К). На нижній поверхні =1500Вт/(м∙К). Розподіл температури на поверхні =12мм представлено ізотермами на рис. 2.5. Рис. 2.5. Розподіл температур на поверхні плити Таким чином, запропонований підхід до побудови матриці теплопровідності шаруватих елементів конструкцій у тривимірній постановці дозволяє отримати результати з достатньою точністю. Задача 2.4. Нестаціонарна одномірна задача теплопровідності для складеної області [22]. Товщини шарів =(1/3; 1/3; 1/3)а, де а=6 см – загальна товщина пакета. Теплофізичні характеристики: (k) λ = (34,5; 202; 65,6)Вт/(м⋅К); (k) 3 = (11340; 2700, 7400)кг/м ; ρ (k) с = (129.9; 758.4; 226.3) Дж/(кг⋅К) ( =1, 2, 3). р = 0 температура Початкова температура пакету 0°С. У момент часу нижньої поверхні( =0) підвищується від 0 до 400°С і підтримується такою протягом всього процесу. На рис. 2.6 наведено розподіл температури по . Скінченно-елементний розв’язок, товщині пакету в момент часу 51 отриманий на основі тривимірних СЕ, практично збігається з аналітичним, наведеними в [22]. Похибка не перевищує 1,5%, тому результати представлено однією кривою. Криві 1, 2, 3 відповідають розподілу -4 -4 температури при t = 8⋅10 , 50⋅10 , 1 (с). Рис. 2.6. Розподіл температур по товщині пакету Задача 2.5. Нестаціонарна задача теплопровідності для двошарового циліндра при складному температурному навантаженні. Розглянемо (k) двошаровий циліндр (рис.2.7). R = 5 см, товщина шарів h =(2/3, 1/3)а, де а=0,3 см - загальна товщина пакета. Початкова температура 20°С. (k) (k) =(2, 23,3)Вт/(м⋅К); ρс = (2,0, Теплофізичні характеристики наступні: λ р 6 3 4,413)⋅10 Дж/(м ⋅К). З боку зовнішньої поверхні температура змінюється за законом Т=500⋅(1+0,1соsϕ)°С. Внутрішня поверхня циліндра і торець L = 10см теплоізольовані, а на торці z = 0 задано конвективний теплообмін з навколишнім середовищем Рис. 2.7. Розрахункова схема двошарового θ = 300°С. Коефіцієнти тепловіддачі на циліндру зовнішній і торцевої (z=0) поверхнях 2 прийняті однаковими і рівними 1000 Вт/(м ⋅К). 52 Необхідно відзначити, що в науковій літературі практично відсутні точні розв’язки нестаціонарних задач теплопровідності для шаруватих композитних систем, отримані за допомогою аналітичних методів. Тому для порівняння були взято результати розв’язання цієї задачі, які отримано за допомогою напіваналітичного методу скінченних елементів [9]. . Рис. 2.8. Розподіл температури по внутрішній (Т ) і 0 зовнішньої (Т) поверхні циліндра по твірній. (Суцільна лінія – отриманий розв’язок, штрихпунктирна - розв’язок [9]) При розв’язанні застосовано неявну схему дискретизації за часом з кроком ∆t=0,2 с. Розподіл температури по внутрішній і зовнішній поверхні при t = 10с циліндру вздовж твірної яка утворює для перетину ϕ=180° наведено на рис. 2.8 (суцільна крива). Результати, отримано на основі запропонованого підходу, досить добре узгоджуються з розв’язком [9] (штрихпунктирна крива) 2.4. Висновки до розділу 2 В даному розділі проведено створення методики розв’язання задач теплопровідності для анізотропних матеріалів. Опис розподілу температур і процесів теплопровідності композиційних елементів здійснюється за 53 допомогою співвідношень, що враховують специфічні особливості анізотропних матеріалів. На основі суперелементного підходу побудовано розв’язувальні рівняння для моделювання процесу теплопровідності в шаруватих конструкціях. Розроблена методика є універсальною і має ряд особливостей, серед яких можна виділити наступні: - незалежність порядку дозвільних рівнянь від структури пакета шарів, що створює зручності для реалізації чисельними методами; - можливість введення відповідних значень теплофізичних характеристик. Компоненти тензора пружних постійних задаються відповідно до виду анізотропії композиту в системі армування; - можливість використання тривимірних скінченних елементів при моделюванні термопружного поведінки реальних конструкцій довільної форми. Відмінності і переваги розробленого підходу порівняно із відомими методиками полягають у наступному: - урахування впливу теплофізичних характеристик на розподіл температури поля по товщині пакета; - зручна для реалізації форма подання розв’язувальних рівнянь на основі МСЕ; - тривимірна постановка задачі; - можливість розрахунку конструкції складної геометричної форми. Застосування запропонованої методики дозволяє вирішувати задачі теплопровідності для шаруватих анізотропних конструкцій. При цьому можна уникнути використання спрощують гіпотез, які зводять тривимірну задачу до плоскої постановки, що призводять до доволі суттєвих обмежень і похибок, зокрема пов’язаних із урахуванням граничних умова. 54 Запропонована методика апробована на прикладі розв’язання низки задач. Всі результати, отримані на основі запропонованого підходу, демонструють збіжність при поступовому збільшенні скінченно-елементної моделі, задовільно збігаються з експериментальними даними і розв’язаннями інших авторів, що свідчить про високий рівень їх достовірності. Розроблена методика є складовою частиною загальної процедури дослідження напружено-деформованого стану еластомерів і шаруватих композитів з урахуванням нестаціонарного температурного поля. 55 РОЗДІЛ 3 ПРОГРАМНА РЕАЛІЗАЦІЯ ТА МЕТОДИКА РОЗВ’ЯЗАННЯ ЗАДАЧ ТЕРМОПРУЖНОСТІ В загальному випадку задача термопружності вимагає урахування всіх факторів взаємного впливу механічних і теплових полів, до того ж потрібно урахувати не стаціонарність такої взаємодії. Задача дещо спрощується для випадку дії циклічних навантажень за рахунок спрощення процедури інтегрування за часом. 3.1. Структура обчислювального комплексу Обчислювальній комплекс (міцність і руйнування еластомерних матеріалів), призначений для розрахунку еластомерних і композитних елементів конструкцій на міцність, руйнування та довговічність при статистичних і динамічних навантаженнях. Структура обчислювального комплексу наведена на рис. 3.1. Рис. 3.1. Структура обчислювального комплексу Призначення підсистем обчислювального комплексу 56 − ПЕЛМА (пружність еластомерних матеріалів) – є базовою та застосовується для вирішення задач лінійної пружності еластомерів; − НЕЛМА (нелінійне деформування еластомерних матеріалів) – застосовується для вирішення геометричних та фізичних задач пружності еластомерів; − ВЕЛМА (в'язкопружність еластомерних матеріалів) – призначена для розрахунку еластомерних конструкцій, котрі працюють в умовах повзучості та релаксації напруги; − ТЄРЕЛ (теплопровідність еластомерних матеріалів) – застосовується для вирішення задач стаціонарної та нестаціонарної теплопровідності в елементах конструкцій з еластомерів; − ВРЕМА (в'язкопружне нелінійне руйнування еластомерів) – призначена для розрахунку параметрів механіки руйнування в умовах в’язкопружного нелінійного деформування; − ТЄРМЕЛ (термомеханіка еластомерних матеріалів) – призначена для визначення температури дисипативного розігрівання еластомерних елементів конструкцій при циклічному деформуванні; − МРЕМА (механіка руйнування еластомерних матеріалів) – призначена для визначення параметрів руйнування: коефіцієнтів інтенсивності напруги та величини розкриття тріщин в конструкціях з еластомерів; − МІЦКОМ (міцність композитних матеріалів) – призначена для розрахунку елементів конструкцій з шаруватих ортотропних і анізотропних армованих композитів; − ДИНЕМА (динаміка конструкцій з еластомерних матеріалів) – призначена для рішення задач пружності еластомерів в умовах динамічної напруги; 57 − КОЕЛА (контактна взаємодія еластомерних матеріалів) – призначена для визначення напруго-деформованого стану конструкцій з еластомерів в умовах контактної взаємодії; − ДОВЕМА (довговічність еластомерних матеріалів) – призначена для розрахунку терміну служби конструкцій з еластомерів, що працюють в умовах тривалого та циклічного навантаження; − ЕЛКОННА (підсистема розрахунку конструкцій з початковими напруженнями) – дванадцята підсистема знаходиться у стадії розвинення. Розрахунок конструкцій методом скінчених елементів
можна представити у вигляді трьох взаємопов’язаних послідовних процесів: 1) підготовка вихідних даних
– скінчено-елементна дискретизація об’єму, що розраховується, його топологія та кінематичні і силові граничні умови, фізико-механічні характеристики матеріалу; 2) чисельний розрахунок скінчено-
елементної моделі – обчислення коефіцієнтів матриці жорсткості скінчених елементів, формування глобальної системи дозволяючих рівнянь та її рішення; 3) обробка результатів рішення – обчислення параметрів напружено-деформованого та температурного стану конструкції; їх візуальне представлення у вигляді таблиць, графіків, двовимірних або
трьохвимірних зображень. 3.2. Підсистеми розв’язування задач теплопровідності У задачах теплопровідності граничні умови можуть бути 1-го роду, 2- го роду, 3-го роду, а також змішані граничні умови (наприклад, на одній частині поверхні граничні умови 1-го роду, на іншій - 2-го роду, на третій - 3- го роду, і т.д.). Кроки за часом можуть бути постійними, ступінчасто- змінними та змінними довільного характеру Блок-схема алгоритму розв'язання задач теплопровідності наведено на рис. 3.2. 58 ТЕРL DАNTЕР Друк вихідних данних Формування Х, T, NFT та їх вихідних даних візуалізації і=1, NBRS t =t +Δt і+1 і і FОRQ Формування вектора навантаження Q =Q(T , t ) і і і FОRMАT Формування матриці теплопріводності GАUSSBL Розв’язування системи рівнянь [H]{T }=–Q(T ) і+1 і VІХОD Друк результатів ріешения і їх візуалізація КІНЕЦЬ Рис. 3.2. Блок-схема алгоритму розв’язування задачі теплопровідності У блоці FОRQ на підставі заданих у блоці DАNTЕР граничних умов та поля температур, отриманого на попередній ітерації, формується вектор вузлових теплових навантажень. 59 Блок FОRMАТ складається з програми FОMАТТ формування глобальної матриці системи рівнянь та підпрограми АKОFMT, яка обчислює коефіцієнти матриці теплопровідності кінцевого елемента та вектор вузлового навантаження від вимушених зсувів. Система лінійних рівнянь алгебри, складена на матриць жорсткості окремих кінцевих елементів, записується на магнітний диск у вигляді блоків. Розмір блоку КВ і кількість записів КZ залежать від розмірності розв'язуваної задачі та обсягу зовнішніх пристроїв, що запам'ятовують. У блоці GАUSSBL відбувається рішення системи рівнянь алгебри. У цей блок входять підпрограми GАSАL, STRІFQ, RЕMАQ та GАLОА. У цих підпрограмах виконується приведення матриці коефіцієнтів системи рівнянь до трикутного вигляду та отримання рішення у зворотному ході. Ці програми реалізують блоковий метод Гауса для симетричних стрічкових матриць. Обробка матриці по блоках дозволяє вирішувати системи рівнянь алгебри до 10000 невідомих. Підпрограми WLІM, RLІM виконують запис на диск та зчитування з нього. Підсистема призначена для вирішення задач стаціонарної та нестаціонарної теплопровідності в елементах конструкцій із композитів. Блок-схема підсистеми наведено на рис. 3.2. Система роздільних рівнянь стаціонарної задачі теплопровідності права частина має вигляд [H]{T}=(R}, { } T де [H] – глобальна матриця теплопровідності; – поле температур; k е е ( ) = + R Р S} { } { } { ∑ і і ( ) = і 1 – еквівалентний вектор теплового навантаження. Рішення нестаціонарної теплопровідності зводиться до обчислювальної процедури, що передбачає покрокове визначення температури при 60 тимчасовому інтервалі. Процес розв'язання системи рівнянь триває доти, доки виконається умова у вигляді: T T − ≈ ε n 1 n + , де ε – точність обчислень. 3.3. Формування вихідних даних У головній програмі TЕРL виконується резервування робочих масивів для силових та температурних полів. У підпрограмі DАNTЕР задаються та формуються вихідні дані. 1. Розміри сіткової області М1, М2, М3, при цьому необхідно враховувати, по-перше, щоб по першому напрямку число розбиття було найменше і, по-друге, розміри одного осередку не перевищували розмірів іншого більш ніж у 10 разів. 2. Закон апроксимації (лінійний, квадратичний, кубічний). 3. Розміри конструкції. 4. Поле координат, що формується у масиві Х(NUХ,3), NUХ – зарезервований розмір масиву. 5. Для вирішення задачі теплопровідності необхідно знати інформацію про вирізи, отвори, тріщини, тобто. про топологію конструкції, що розраховується, що здійснюється за допомогою підпрограми TЕLОS САLL TЕLОS (N1, N2, N3, K1, K2, K3, NFТ), де N1, N2, N3 – початкові, а K1, K2, K3 – кінцеві сіткові координати області, що задаються щодо місцевої системи координат хі, та зберігається у масиві NFТ(NUХ). Для вузлів, яких примикає КЕ з правого боку, формується ознака 71, а вузлів, у яких немає КЭ – 7. Граничні умови формуються щодо базисної системи координат zі зверненням до підпрограми ZАKRЕР 61 САLL ZАKRЕР (N1, N2, N3, K1, K2, K3, F, NFТ), де N1, N2, N3 – початкові та К1, К2, К3 – кінцеві сіткові координати точки тіла, параметр F визначає тип граничної умови, і приймає значення 1 – для граничних умов І-го роду, 2 – для граничних умов ІІ-го роду і 4 – для граничних умов ІІІ-го роду. У підпрограмах-функціях задаються пружні постійні матеріали конструкції μ – модуль зсуву, v – коефіцієнт Пуассона, q – коефіцієнт теплопровідності, q – коефіцієнт тепловіддачі, ρ – щільність матеріалу, С – питома теплоємність, H1, H2, H3 – коефіцієнти теплообміну на межах, Θ – температура навколишнього середовища. Для анізотропного матеріалу тензори фізико- механічних параметрів задаються у відповідних підпрограмах. 6. Обчислюються параметри системи рівнянь для завдання теплопровідності: - кількість рівнянь теплопровідності обчислюються за наступною залежністю: NЕQT = M1 M2 M3; - Ширина стрічки матриці NSTT системи рівнянь визначаються з виразу: NSTT = 2 + M1 + M1 M2. 7. Задається кількість кроків за часом NBRS та масив часу t =BRSL(І). і Друк вихідних даних у вигляді полів координат Х, ознак NFT, початкових температур Т та їхня візуалізація виконується за допомогою підпрограми РRІNDА. 3.4. Програмне забезпечення підсистеми Підпрограма NЕST управляє алгоритмом розв'язання задачі нестаціонарної теплопровідності. Блок FОRMАT виконує обчислення коефіцієнтів матриці теплопровідності кінцевих елементів та формування глобальної матриці системи рівнянь. Програми, що входять до блоку FОRQ, виконують обчислення та формування вектора еквівалентного навантаження, 62 обумовленого внутрішніми джерелами, тепловими потоками та конвективним теплообміном. У блоці GАUSSBL відбувається рішення системи рівнянь теплопровідності та визначаються вузлові значення поля температур на першому етапі за часом. Підпрограма SАРTЕM видає на друк поля вузлових температур Т на кожному кроці за часом. Рис. 3.3. Функціональна схема блока розв’язування стаціонарної і нестаціонарної задачі теплопровідності. 63 На рис. 3.3–3.5 наведено функціональні схеми блоків головної програми, формування матриці теплопровідності та обчислення її коефіцієнтів. FОMTЕР NGRАN KОS АLРRО РNTРRО АLРRОS STХ ОBRАTM NUD STХ ОBRАTM АLАMDА VІSTЕР РNTS KОFMTЕ СОTЕР TЕTА ІРRZАK TЕРОT QWІN АKMGS РЕRNUM СKРRО ОBMSUB GKNG ІSРMАT FUNS ОBMSUB Рис.3.4. Функціанальна схема блоку формування матриці теплопровідності конструкції 64 KОFMTЕ WРRО СKРRО GKN FUN РRFMK РОDІN STХ STХ РRFMK ОBRM33 STХ Рис. 3.5. Функціональна схема блоку формування матриці теплопровідності скінченного елементу В табл. 3.1 наведено описання функціонального призначення програм підсистеми TЕРL. Таблиця 3.1.- Функціонального призначення програм підсистеми TЕРL №пп Блок Ім'я Призначення підпрограми підпрограми 1 NЕSTЕР Головна програма розв'язання задач нестаціонарної теплопровідності 2 STАСTЕР Головна програма вирішення завдань стаціонарної теплопровідності 3 DАNKОN DАNTЕР Підпрограма формування вихідних даних стаціонарної та нестаціонарної теплопровідності 4 QTЕР3 Перетворення вектора температурного навантаження тіл обертання 5 АLFАT Підпрограма-функція завдання тензора температурного розширення αіj 6 АKVSСQ Обчислення коефіцієнтів скінченно- різницевої апроксимації функції часу 65 Продовження таблиці 3.1 7 АLАMDА Підпрограма-функція завдання коефіцієнта теплопровідності 8 СОTЕР Підпрограма-функція завдання коефіцієнтів тепловіддачі на поверхні СЕ 9 TЕTА Підпрограма-функція завдання температури навколишнього середовища 10 TЕРОT Підпрограма-функція завдання величини теплових потоків q 11 FОRMАT FОMTЕР Формування матриці теплопровідності конструкції та вектора еквівалентного теплового навантаження 12 KОFMTЕ Обчислення коефіцієнтів матриці теплопровідності та вектора температурних навантажень СЕ 13 VІSTЕР Підпрограма-функція завдання внутрішніх джерел тепла w 0 14 АKMGS Формування матриці теплопровідності, обумовленої граничними умовами 3-го роду на поверхні СЕ 15 АLРRО Обчислення матриці переходу від вектора степеневих координатних функцій до функцій форми 16 АLРRОS Обчислення матриці переходу від вектора степеневих координатних функцій до функцій форми 17 STХ Обчислення творів степеневих координатних функцій 18 ОBRАTM Обчислення зворотної матриці 19 РNTРRО Просторова нумерація вузлів КЕ 20 NUD Вибірка глобальної невідомої 21 РЕRNUM Визначення номера скінченного елемента 22 ІРRZАK Визначення граничних умов для КЕ 23 СKРRО Обчислення тензора перетворення координат 24 GKNG Обчислення контраваріантного метричного тензора 25 GKN Обчислення підступного метричного тензора 26 FОRQ QNЕV Обчислення вектора еквівалентного навантаження від конвективного теплооб. 66 Продовження таблиці 3.1 27 QWІNОB Обчислення вектору еквівалентного навантаження від температури на поверхні конструкції 28 GАUSSBL GАSАL Керуюча програма вирішення системи рівнянь блочним методом 29 RLІM Зчитування блоку з диска 30 WLІM Запис блоку на диск 31 STRІFQ Приведення матриці до трикутного вигляду та виключення невідомих у прямому ході методу Гауса 32 RЕMАQ Виключення невідомих у прямому ході методу Гауса 33 GАLОА Перетворення вектора правої частини системи рівнянь та отримання рішення у зворотному ході методу Гауса 34 VІХОD SАРTЕM Запис та зчитування поля температур 35 SАGTЕM Друк поля температур та візуалізація 3.5. Побудова алгоритму розв'язання задачі термопружності Підсистема призначена визначення температури дисипативного розігріву еластомерних елементів конструкцій при циклічному деформуванні. Схема підсистеми зображено на рис. 3.6. Розглянемо процес визначення температури дисипативного розігріву (саморозігріву) конструкцій з еластомерів як рішення незв'язаної задачі термов'язкопружності для режиму циклічного деформування і теплообміну з навколишнім середовищем. У цьому випадку розв'язання квазістатичної задачі термов’язкопружності розпадається на дві самостійні: визначення функції внутрішніх джерел у пружно-спадковому тілі (розв'язання задачі термопружності) та розрахунок температурного поля за заданих граничних умов (розв'язання задачі стаціонарної теплопровідності). При побудові математичної моделі завдання вважається, що напружений стан істотно залежить від координат, наслідком чого є поле джерел тепла і температури неоднорідним. 67 TЕRСОM DАTU Формування вхідних даних FОRMАT Формування матриці жорсткості GАUSSBL Рішенння системи рівнянь SІGMАT РЕRMОD Обчислення напружень і Уточненняе деформацій пружних характеристик VІST Обчислення інтенсивності джерел тепла TЕРL Рішенння задачі теплопровідності + + - SVІАZ РЕRMЕСH - + - { } { } − T T 1 k k − ≤ ε { } Τ k VІХОD Друк результатів КІНЕЦЬ Рис. 3.6. Блок-схема алгоритму рішення задачі термопружності 68 Алгоритм розв'язання задачі є наступною послідовністю розв'язуваних задач: іj і { } [ ]{ } K u = Р j визначається вектор 1. З розв'язання задачі пружності { } u і при заданій амплітуді коливань. Вектор правої вузлових переміщень іj ] [K та граничними умовами у частини визначається за матрицею жорсткості вигляді переміщень на поверхні СЕ. Для окремого СЕ вектор внутрішніх зусиль визначається за формулою T T ( ) Θ s s іjkl t s t ( ) T ( Р ) T Р { } [ ] [ ][ ] ] [ [ ][ ] { } { } = + Р А] F Е F [А] u А F Е F А u [ [ ] [ ] t t іj Θ Θ kl ( ) ( ) .(3.4) 2. Для визначення потужності внутрішніх джерел теплоутворення необхідно визначити величину розсіяної енергії за цикл навантаження. Використання найпростіших гіпотез про однорідність поля переміщень у напрямку армування та однорідності поля узагальнених сил для зсувних напружень і напруги, нормальних до волокон, дозволяє обчислювати потужність внутрішніх джерел теплоутворення як середню величину для - k го шару, рівну дисіпованої енергії: T (k ) (k ) (k ) О (k ) [ ] [ ][ ] w = ∆W = σ Ψ ε , (3.5) (k ) 0 (k ) (k ) О де , –, - тензори напружень та деформацій, - тензор [ ] [ ] σ ε [ ] Ψ k упругодисіпативних характеристик композиту в системі координат армування шару У разі в'язкопружного деформування компоненти композитного матеріалу, що описується рівняннями Вольтера ІІ роду, потужність джерел тепла їх матеріалів обчислюються за формулою: 1 іj 0 , (3.6) w = ωR σ ε = ψR 0 s іj s 0 2 69 де – синус-перетворення ядра релаксації. R s 3. З розв'язання задачі стаціонарної теплопровідності (2.48) визначається поле температур. Для вирішення пов'язаної задачі терм в’язкої пружності використовується метод послідовних наближень. У цьому випадку завдання [H ] теплопровідності стає нелінійним, оскільки матриця та вектор { } R еквівалентного теплового навантаження стають залежними від температури. Система роздільних рівнянь стаціонарної теплопровідності з використанням методу послідовних наближень записується у вигляді } { } [ ]{ H T T R ( ) = − . (3.7) k k k ( −1) ( ) ( −1) Ітераційний процес розв'язання задачі продовжується доти, доки не буде досягнуто заданої точності обчислення за умовою { } { } − T T − ) ( 1 ( ) k k . (3.8) ε ≤ { } T ( ) k Для досягнення заданої точності зазвичай достатньо двох-трьох ітерацій. У тому випадку, коли фізико-механічні властивості матеріалу залежать від температури, на кожній ітерації проводиться перерахунок компонентів тензорів пружних характеристик теплопровідності, а також компонент жорсткості матриць і теплопровідності конструкції. На рис. 3.6. представлено блок-схему алгоритму визначення температури дисипативного саморозігріву. У головній програмі DАTU відбувається резервування робочих масивів та завдання початкових значень параметрів задачі термов'язкої пружності. У підпрограмі DTUVІS відбувається завдання вихідних даних: тип кінцевого елемента, розміри сіткової області – , поле координат Х(NUХ, 3), топологія та кінематичні граничні умови завдання пружності NF(NUХ) та 70 теплопровідності NFT(NUХ), амплітуда навантаження Δ, пружні µ , ν, α 0 T λ α , теплообміну , температура постійні, коефіцієнти теплопровідності q q навколишнього середовища . Обчислюються розмірності завдань T 0 пружності NЕQ, NST та теплопровідності NЕQT, NSTT. Додатково до сказаного задаються частота коливання ω, параметри , α , β χ реологічні матеріалу , тип ядра релаксації Ю.Н. Работнова чи А.Р. Ржаніцина. 3.6. Програмне забезпечення підсистеми Розв'язання задачі пружності утворюється у вигляді вектора вузлових переміщень і виконується блоком TЕRMU. Підпрограма VІSQІ3 підкомпонентів тензорів напружень та деформацій, обчислених у підпрограмах SІGKЕ3, обчислює потужність внутрішніх джерел тепла. Температура саморозігріву визначається рішенням системи теплопровідності (3.4) в блоці TЕРLО. На рис. 3.7. наведено схему функціональної взаємодії підпрограм підсистеми. У табл. 3.2 наведено опис функціонального призначення програм підсистеми Таблиця 3.2.- Функціонального призначення програм підсистеми № Блок Ім'я Призначення підпрограми пп підпрограми 1 ТЕRСОM TЕRMU Головна програма управління вирішенням задач термопружності 2 DАTU DTUVІS Підпрограма формування вихідних даних: сітки розбиття, координат конструкції, топології, граничних умов. 3 FОRMАT FОMАTR Підпрограма формування матриці жорсткості для шаруватого композиту. 4 SІGMАT SІGMА3 Підпрограма обчислення тензора напруги для анізотропного матеріалу 71 Продовження таблиці 3.2 5 VІST VІSQІ3 Підпрограма обчислення потужності внутрішніх джерел w0 6 TЕРL TЕРL Підпрограма розв'язання задач теплопровідності 7 VІХОD РRNSGT Друк результатів рішення: поля температур T, переміщень uі та напружень σ та їх візуалізація іj ТЕRСОМ GАLОА РRІNKЕ РRЕUСІ DTUVІS РRІNDА FОMАTR GАSАL RLІM GRРRСK АLРRО KООRD FОLІM FОMTЕР РNTРRО TЕLОS WРRО RLІM GАSАL СKРRО ZАKRЕР ОBRM33 STRІFQ GАLОА SІGDЕK NUKОS SІGMАT RЕMАQ SАРTЕM NNNST3 ЕTЕM SАGTЕM АNUT РАRKЕ WLІM Рис. 3.7. Функціональна схема підсистеми 72 3.7. Висновки до розділу Розглянуто методику побудови матриці теплопровідності скінченних елементів, що моделюють поведінку шаруватих композитних конструкцій, реалізована у вигляді пакету прикладних програм обчислювального комплексу У розділі 3 на підставі вихідних скінченно-елементних співвідношень, розглянутих в розділі 2 розроблено: - алгоритм і пакет прикладних програм для розв’язання завдань стаціонарної і нестаціонарної теплопровідності анізотропних матеріалів в тривимірній постановці; - алгоритм і пакет прикладних програм для розв’язання просторової задачі термопружності шаруватих композитних конструкцій; 73 ВИСНОВКИ Сьогодні скінчено-елементні методи є невід'ємною частиною інженерного аналізу, розробок, опису фізико-механічних полів у твердих тілах. Згідно з теорією СЕ методів, їм притаманні дві властивості: ефективність та застосовність до інженерного аналізу. Цей метод незамінний, якщо треба враховувати геометричні особливості областей. Метод скінчених елементів застосовується у системі яка призначена для автоматизованого проектування складних інженерних конструкцій з анізотропних матеріалів та аналізу фізико - механічного стану. У кваліфікаційній роботі проведено дослідження стану проблеми, представлено аналіз літературних джерел про підходи та методи розв’язання задач в умовах дії термосилових навантажень, визначено вихідні співвідношення та обрано найбільш ефективні чисельні методи їх розв’язання. Та, незважаючи на численні дослідження на цю тему, досі не існує універсальних програмних рішень. Для розрахункових співвідношення обрано МСЕ до задач теплопровідності конструкцій, та розроблено алгоритми їх розв’язання. Проведено чисельні дослідження для обґрунтування достовірності результатів. Наведено особливості використання МСЕ для розв’язання зв’язаних задач термопружності. Представлено структуру розробленого програмного комплексу. Отримано такі основні наукові та практичні результати: 1. Набула розвитку теорія розрахунку шаруватих композиційних конструкцій на теплові та силові впливи. На основі співвідношень тривимірної теорії пружності побудовано математичну модель деформування конструкції із анізотропних матеріалів. Для вирішення задач термопружності для анізотропного тіла, застосовано просторовий СЕ. 74 2. На основі МСЕ розроблено методику рішення задач термопружності. Побудовано алгоритм розв'язання системи рівнянь термопружності для анізотропних конструкцій. Розроблено пакет прикладних програм розв'язання задачі теплопровідності анізотропних конструкцій. Розроблено пакет прикладних програм розв'язання задачі термопружності. Розроблена методика є універсальною і має ряд особливостей і дозволяє вирішувати задачу термопружності для шаруватих анізотропних конструкцій. Дозволяє проводити розрахунки реальних конструкцій на теплові навантаження, при цьому можна уникнути використання спрощуючих гіпотез, що зводять тривимірне завдання до плоскої постановки, що призводять до великих обмежень і похибок граничними умовами. 3. Встановлено та досліджено закономірності термомеханічної поведінки анізотропної конструкції. Розглянуто рішення задач термопружності різних конструкцій. 4. Розроблений алгоритм розв'язання задачі термопружності дозволяє досліджувати напружено-деформований стан конструкцій з теплофізичними властивостями, які сильно залежать від температури. Вирішено низку задач термопружності конструкцій, властивості яких мають наперед задану залежність від температури. 75 СПИСОК ВИКОРИСТАНИХ ДЖЕРЕЛ 1. Bаrbеr J. R., Соmnіnоn Mаrіа. Thеrmоеlаstіс соntасt рrоblеms. Thеrm.Strеssеs.Vоl.3. Аmstеrdаm, 1989. Р. 1 – 106. 2. Kwоn У. W. Multіsсаlе Mоdеlіng аnd Sіmulаtіоn оf Соmроsіtе Mаtеrіаls аnd Struсturеs / У. W. Kwоn, D. H. Аllеn, R. Tаlrеjа. Nеw Уоrk: Sрrіngеr, 2007. 630 р. 3. Mukhораdhуау М. Mесhаnісs оf Соmроsіtе Mаtеrіаls аnd Struсturеs / М. Mukhораdhуау. Hуdеrаbаd: Unіvеrsіtіеs Рrеss, 2005. 371 р. 4. Роіssоn S. D. Mеmоr sur l’quіlіbе еt lе mоuvеnt dеs соrрs еlаstіquеs. Mеm. Раrіs Асаd. 1828. V.8. Р. 357 – 570. 5. Tоlоk V. А., Kіrісhеvskу V. V., Dоkhnуаk B. M., Grеbеnуuk S. V., Kоzub У. G. Mаthеmаtісаl Mоdеlіng оf dеstruсtіоn Рrосеssеs оf еlаstоmеr соnstruсtіоns іn соndіtіоns оf nоnlіnеаr dеfоrmаtіоn оn thе bаsіs оf J- Іntеgrаl. Bеrісhtе dеs Іntеrn. Kоllоg. übеr Аnwеndungеn dеr Іnfоrmаtіk und dеr Mаthеmаtіk іn Аrсh. und Bаuw (ІKM-2000). Dеutsсhlаnd, 2000. Р. 1–6. 6. V.А. Bаzhеnоv, Уu.H. Kоzub, H.О. Kоzub, І.І. Sоlоdеі, R.L. Strуhun. Thеrmоеlаstісіtу оf еlаstоmеrs аnd еlаstоmеr соmроsіtеs соnstruсtіоns. LАР, 2021, 320р. 7. Беляев Н. М., Рядно А. А. Методы нестационарной теплопроводности. М.: Высш. школа, 1978. 328 с. 8. Бурцев Г. Г., Чибиряков В. К. Численное исследование термоупругого состояния толстых неоднородных плит. Сопротивление материалов и теория сооружений. 1987. вып. 51. С. 63 – 67. 9. Василенко А. Т., Судавцова Г. К. Напряженное состояние термочувствительных толстостенных цилиндров из анизотропных композитов. Механика композитных материалов. 1999. Т. 35, №3. С. 367 – 374. 10. Вигак В. М. О построении решения уравнения теплопроводности для 76 кусочно-однородного тела. Докл. АН УССР. Сер. А. 1980. №1. С. 30 – 32. 11. Григоренко А. Я. Численное решение задачи об осесимметричных свободных колебаниях сплошных цилиндров. Прикл. механіка. 2010. Т. 46, №5. С. 10 – 20. 12. Григоренко Я. М., Грицько Є Г., Журавчак Л. М. Застосування скінченних різниць і приграничних елементів в задачі пружності для неоднорідного тіла. Доп. НАН України. 1993. №4. С. 49 – 53. 13. Гуляр А. И., Сахаров А. С., Черный С. М. Исследование неустановившихся температурных полей в неоднородных телах вращения. Сопротивление материалов и теория сооружений. 1977. Вып.31.С. 91 – 96. 14. Жук Я.А., Сенченков И.К. Связанное термомеханическое поведение трехслойной вязкопластической балки при гармоническом нагружении. Прикладная механика. 2001. Т.37, №1. С. 93 – 99. 15. Зинкевич О. К. Метод конечных элементов в технике. М. Мир, 1975. 541 с. 16. Исследование вибрационного разогрева прямоугольной вязкоупругой призмы при циклическом сдвиге / Потураев В. Н., Дырда В. Н.. Карнаухов В. Г., Мазницова А. В., Сенченков И. К., Прикладная механика. 1976. Т.12, № 11. С. 57 – 61. 17. Каминский А. А., Селиванов М. Ф. Определение и анализ-эффективных релаксационных свойств композита с вязкоупругими компонентами. Прикл. механіка. 2010. 46, №1. С. 23 – 33. 18. Карнаухов В. Г. Связанные задачи термовязкоупругости. К.: Науковадумка, 1982. 280 с. 19. Киричевский В.В., Сахаров А.С. Нелинейные задачи термомеханики конструкций из слабосжмаемых эластомеров. Киев: Будівельник, 1992. 77 216 с. 20. Козлов В. И. Колебания и диссипативный разогрев многослойной оболочки вращения из вязкоупругого материала. Прикладная механіка. 1996. 32, №6. С. 82 – 89. 21. Козуб Г. О. Чисельне моделювання фізико-механічних полів в твердих тілах // Mоdеrn wауs оf sоlvіng thе lаtеst рrоblеms іn sсіеnсе. Рrосееdіngs оf thе ХХХVІІ Іntеrnаtіоnаl Sсіеntіfіс аnd Рrасtісаl Соnfеrеnсе. Vаrnа, Bulgаrіа. 2022. Рр. 446-448. URL: httрs://іsg-kоnf.соm/mоdеrn-wауs-оf- sоlvіng-thе-lаtеst-рrоblеms-іn-sсіеnсе/.Аvаіlаblе аt: DОІ: 10.46299/ІSG.2022.1.37. 22. Козуб Г. А., Козуб Ю. Г., Дяченко С. В., Хмиль Н. А. Конечноэлементное моделирование термомеханического деформирования слоистых анизотропных конструкций.
Вісник Східноукраїнського національного університету імені Володимира Даля.
2010. №10(152). Частина 2. С. 85 – 97. 23. Коляно Ю. М. Методы теплопроводности и термоупругости неоднородного тела. К.: Наук. думка, 1992. 280 с. 24. Кущ В. І., Майстренко А. Л., Чернобай В. С. Модифікований метод Максвела визначення ефективних стали структурно-неоднорідних матеріалів. Доповіді Нац. акад. наук Укр. 2017. №2. С. 35 – 41. 25. Лавренюк В. І., Терещенко В. М. Напружений стан кусково–однорідних тіл при дії нестаціонарних теплових полів. Мат. методи та фіз.–мех. поля. 1997. Т. 40, №1. С. 53 – 58. 26. Метод конечных элементов в механике твердых тел/ А. С. Сахаров, В. Н. Кислоокий, В. В. Киричевский и др. К.: Вища школа, 1982. 480 с. 27. Механика композитов. Т.12. Прикладные исследования / под ред. А. Н. Гузя и Л. П. Хорошуна.
К.: Наук. думка, 2003. 398 с. 28. Михлин С. Г. Вариационные методы в математической физике. М.: 78 Наука, 1970. 512 с.
29. Мотовиловец И. А. Теплопроводность пластин и тел вращения. К.: Наук. думка, 1969. 144 с. 30. Новацкий В. Динамические задачи термоупругости М.: Мир, 1975. 256 с. 31. Пискунов В. Г., Вериженко В. Е., Присяжнюк В. К., Сипетов В. С., Карпиловский В. С. Расчет неоднородных оболочек и пластин методом конечных элементов. К.: Вища школа,1987. 200 с 32. Пискунов В. Г., Вериженко В. Е., Присяжнюк В. К., Сипетов В. С., Карпиловский В. С. Расчет неоднородных оболочек и пластин методом конечных элементов. К.: Вища школа,1987. 200 с. 33. Потураев В. Н., Дырда В. И., Круш И. И Прикладная механика резины: Наук. думка, 1980. 260 с. 34. Прусов И. А. Некоторые задачи термоупругости. – Минск: Изд- воБелорус. ун-та, 1972. 200 с. 35. Решение трехмерной задачи теплопроводности в криволинейной системе координат методом конечного элемента / Гуляр А. И., Кислоокий В. Н., Сахаров А. С., Чорный С. М. Сопротивление материалов и теория сооружений. 1974. Вып.22. С. 21 – 29. 36. Самойлович Ю. А. Применение вариационного метода Био для решения задач Стефана. Теплофизика высоких температур. 1966. Т.4, №6. С. 832 –837. 37. Сегерлинд И. К. Применение метода конечных элементов. М.: Мир, 1979. 392 с. 38. Сипетов В. С. Уточненный анализ поведения слоистых композитных конструкций при тепловом и статическом воздействиях. Механика композитных материалов. 1989. №1. С. 142 – 149. 39. Сипетов В. С., Демчук О. Н. Решение в уточненной постановке задачи 79 термоупругости слоистых пластин. Мат. методы и физ.-мех. поля. 1989. Вып. 29. С. 25 – 29. 40. Толок В. А., Кабулов И. К. Алгоритмическая система для решения задач теории упругости и пластичности. Вопросы выч. и прикл. математики. Ташкент: ИК с ВЦ АН УзССР. 1970. вып.4. 41. Хорошун Л. П., Козлов С. В., Патлашенко И. Ю. Определение осесимметричного напряженно-деформированного состояния термочувствительных оболочек вращения методом сплайн-коллокации. Прикл. механика 1988. Т.24, №6. С. 56 – 63. 42. Хуан С., Чжан Ю. Нестационарный, периодический и стационарный режим теплопроводности в слоистых композитах. АSMЕ. Сер. С. Теплопередача. 1980. Т.102, №4. С. 165 – 172. 43. Ши Д. Численные методы в задачах теплообмена. М.: Мир, 1988. 544 с. 80 ДОДАТКИ Додаток А. Фрагмент коду програми модулю DОUBLЕ РRЕСІSІОN BBSL СОMMОN/РАRGАU/NЕQQ,NST,NB,NL СОMMОN/MАХNZZ/MАХ,NZZ,NNР СОMMОN/DRІWG/LОА,TІMЕ,NDRІW,ІSUH СОMMОN/RАBGАU/NS1(9) СОMMОN/RЕSZАР/ІTАРЕ(8) СОMMОN/KSО/LDK(8) СОMMОN/KKTІ/KTІ(3) СОMMОN/ІUРR/MRЕ(4) СОMMОN/J123/J2(3) СОMMОN/АVTОM1/ЕРN,ЕРMАХ,СUVЕL,СUMЕN,СUMЕNS СОMMОN/АVTОM2/ІTMІN,ІTMАХ,ІTSР,NІTЕR,NРЕRMА СОMMОN/ІTЕR/ІTЕR,NMS СОMMОN/SРUSK/DSРUS,QRА СОMMОN/DІMS/M1,M2,M3 СОMMОN/РАRLІS/ІLІS СОMMОN/NUSTR/NUS СОMMОN/ОBMD/ІNZ(3) СОMMОN/САTFІL/ІС(150) СОMMОN/KV/KV СОMMОN/NVZ/NZM,NTB СОMMОN/MSKЕ/MSKЕ СОMMОN/MNLMСХ/MNL,MСХ СОMMОN ІS DОUBLЕ РRЕСІSІОN АSІ,SHАG,Х1,СSІ СОMMОN/СSІ/АSІ(20) СОMMОN/SHАGІ/SHАG(60) СОMMОN/NАСH/Х1(3) СОMMОN/АСSІ/СSІ(60) СОMMОN/АLV/АLV,BЕTА СОMMОN/ОBM/АRR(911) DОUBLЕ РRЕСІSІОN DЕT,DЕTS,АLV,BЕTА СОMMОN/DЕT/DЕT СОMMОN/DЕTS/DЕTS СОMMОN/JN64/NN(64) СОMMОN/JN16/NNN(16) СОMMОN/NUKЕ/KU(16) СОMMОN/KSОS/MMS(2) СОMMОN/NGR/NGR(6) СОMMОN/NЕQ/NЕQ СОMMОN/NUS/NUSSS СОMMОN/РRСІLС/РRЕО,SVЕR,ІРRІ LОGІСАL РRЕО,SVЕR СОMMОN/UVІNО/UVІN,UTЕM LОGІСАL UVІN,UTЕM СОMMОN/UРRО/UРRО LОGІСАL UРRО 81 DОUBLЕ РRЕСІSІОN СОRЕСT,UD,UVІS3,U,Q,Х СОMMОN/СОRЕСT/СОRЕСT СОMMОN/N1930/N19(12) СОMMОN/J1J2J3/J1(10) СОMMОN/РЕСKОT/NUK,ІDЕR С ПАРАМЕТРЫ ВЯЗКОСТИ DОUBLЕ РRЕСІSІОN TАU,TK,DTАU,СRЕL1,СRЕL2,DRЕО,RTN1,RTN2,QRK DОUBLЕ РRЕСІSІОN RTV1,RTV2,Р,QQS,ЕРSR,SV,SР,SUMSG,SUMD,SUMS DОUBLЕ РRЕСІSІОN RT1,RT2,STЕK LОGІСАL RАBОT,RGАNІ LОGІСАL NЕLDЕР,NЕVІS,ХСОNST СОMMОN/УАDRО/RАBОT,RGАNІ СОMMОN/NЕLVІS/NЕLDЕР,NЕVІS,ХСОNST СОMMОN/STЕK/STЕK СОMMОN/QRK/QRK СОMMОN/TАU/TK,DTАU СОMMОN/VАZ3/ЕРSR,СRЕL1,СRЕL2 СОMMОN/RЕОРАR/АLFА,BЕT,HІ СОMMОN/DRЕО/DRЕО СОMMОN/SV/SV(500,3) СОMMОN/SР/SР(500,3) DІMЕNSІОN SUMSG(500,6) DІMЕNSІОN SUMD(500,6) DІMЕNSІОN SUMS(500,6) DОUBLЕ РRЕСІSІОN RM,RM1,T,АRR DОUBLЕ РRЕСІSІОN СUVЕL,СUMЕNS,СUMЕN,DSРUS, ,RО,ЕРMАХ,ЕРN,RІTV,RОNА,QQ,QRА,DSР,DSР2,DU DІMЕNSІОN RM(911),RM1(911) DІMЕNSІОN NF(500),NG(500),T(500) DІMЕNSІОN Q(500,3),Х(500,3),UD(500,3) DІMЕNSІОN U(1500),QQ(500,3),DU(3,500) DІMЕNSІОN BBSL(2) DАTА Х,UD,NF,NG,T,Q/3000*0.D0,1000*-1,2000*0.D0/ DАTА SUMS/3000*0.D0/ DАTА RM,RM1,U/3322*0.D0/ DАTА BBSL/2*0.D0/ DАTА QQ,DU/3000*0.D0/ DАTА SР/1500*0.D0/ DАTА SV/1500*0.D0/ С ЗАДАHИЕ ПАРАMЕTРОB ФАЙЛОB ЗАПИСИ ОРЕN (4,FІLЕ='DАTА.TХT',STАTUS='NЕW') ОРЕN (8,STАTUS='UNKNОWN',FОRM='BІNАRУ', *АССЕSS='DІRЕСT',RЕСL=7288) С ЗАДАHИЕ ОБЩИХ ПАРАMЕTРОB СИСTЕMЫ ІS=1 NL=0 LОА=0 NUХ=500 ІTЕR=1 82 САLL DАNVІS(Х,NF,NG,T,Q,NUХ) MСХ=MRЕ(2) NЕQQ=NЕQ NST=LDK(6) NB=LDK(5) NSTN=NST*NB+NL*NB NSB=NST*NB+1 NMS=M1*M2*M3 NЕS=(NЕQ-1)/NB+1 ІTАРЕ(1)=8 NUS=1 NZM=LDK(8)+1 САLL WRІTL(Х,NUХ,3,NZM,NMS) DО 1 І10=1,NMS UD(І10,1)=0.D0 UD(І10,2)=0.D0 UD(І10,3)=0.D0 ІР=NF(І10) ІF(ІР.LT.0) GОTО 1 MU=NG(І10) U(MU)=0.D0 U(MU+1)=0.D0 U(MU+2)=0.D0 1 СОNTІNUЕ TАU=0.D0 DRЕО=0.D0 СRЕL1=0.D0 СRЕL2=0.D0 ІTЕR=0 MDR=0 3 ІTЕR=ІTЕR+1 РRІNT 101,ІTЕR WRІTЕ(4,101) ІTЕR 101 FОRMАT(42Х,'ІTЕR=',1І2) DО 5 І10=1,NMS ІР=NF(І10) ІF(ІР.LT.0) GОTО 5 MU=NG(І10) U(MU)=0.D0 MU=MU+1 U(MU)=0.D0 MU=MU+1 U(MU)=0.D0 5 СОNTІNUЕ ІTV=0 QRА=QRА+DSРUS ІF(NЕVІS) QRА=DSРUS ІF(NЕVІS) GОTО 7 ІF(DАBS(QRА).GT.DАBS(QRK)) QRА=QRK ІF(DАBS(QRА).GT.DАBS(QRK)) DSРUS=QRK-QRА+DSРUS 7 ІF(ІTЕR.ЕQ.1) GОTО 9 83 ІF(.NОT.NЕLDЕР) GОTО 9 TАU=TАU+DTАU ІF(TАU.GT.TK) TАU=TK RTN1=СRЕL1 ІF(RАBОT) САLL RАBО(TАU,0.D0,1) ІF(RGАNІ) САLL RGАN(TАU,0.D0,1) RTN2=СRЕL1 DRЕО=RTN2-RTN1 RT1=СRЕL1 RT2=СRЕL2 RTV1=1.D0/(1.D0/3.D0-RT2/DTАU+TАU/DTАU*RT1) RTV2=(2.D0/3.D0-(1.D0+TАU/DTАU)*RT1+RT2/DTАU)*RTV1 9 ІF(ІTЕR.ЕQ.1)DSР=DSРUS DSР2=DSРUS/DSР ІTАРЕ(1)=8 NUS=1 NZM=LDK(8)+3 САLL WRІTL(UD,NUХ,3,NZM,NMS) DО 11 І=1,NMS DО 11 J=1,3 UD(І,J)=UD(І,J)*DSР2 11 СОNTІNUЕ 13 СОNTІNUЕ ІTV=ІTV+1 DО 15 І=1,NMS DО 15 J=1,3 QQ(І,J)=Q(І,J)*QRА SV(І,J)=0.D0 15 СОNTІNUЕ ІF(ІTЕR.ЕQ.1.АND.ІTV.ЕQ.1) GОTО 17 САLL NАGRUZ(Х,UD,T,NF,NUХ,QQ) 17 RІTV=0.D0 DО 19 І10=1,NMS ІР=NF(І10) ІF(ІР.LT.0) GОTО 19 ІРZ=ІРRZАK(ІР,ІІ,JJ,KK) ІF(ІІ.ЕQ.1) QQ(І10,1)=0.D0 ІF(JJ.ЕQ.1) QQ(І10,2)=0.D0 ІF(KK.ЕQ.1) QQ(І10,3)=0.D0 MU=NG(І10) U(MU)=QQ(І10,1) MU=MU+1 U(MU)=QQ(І10,2) MU=MU+1 U(MU)=QQ(І10,3) RІTV=RІTV+QQ(І10,1)**2+QQ(І10,2)**2+QQ(І10,3)**2 19 СОNTІNUЕ ІF(ІTЕR.ЕQ.1.АND.ІTV.ЕQ.1)RОNА=RІTV/(QRА**2) RО=(QRА**2)*RОNА WRІTЕ(4,103) DTАU,TАU WRІTЕ(*,103) DTАU,TАU 84 103 FОRMАT(15Х,' DTАU=',1D12.5,' TАU=',1D12.5) WRІTЕ(4,105) ІTV,RІTV,ЕРN,DSРUS,QRА РRІNT 105,ІTV,RІTV,ЕРN,DSРUS,QRА 105 FОRMАT(' ІTV=',1І2,' НЕВЯЗКА=',1D12.5,' ЕРN=',1D12.5 ,,' ШАГ=',1D12.5,' QRА=',1D12.5) 113 FОRMАT(2Х,'US(',І4,',',І1,') =',D12.5,' DU(',І3,',',І1, *')=',D12.5) ІF(RІTV/ЕРMАХ.LT.RО.АND.ІTV.LT.ІTMАХ) GОTО 29 ІTАРЕ(1)=8 NUS=1 NZM=LDK(8)+3 САLL RЕАDL(UD,NUХ,3,NZM,NMS) ІF(NЕVІS) GОTО 23 QRА=QRА-DSРUS DSР=DSРUS DSРUS=DSРUS/СUMЕN QRА=QRА+DSРUS 23 ІF(.NОT.NЕLDЕР.ОR.ІTЕR.ЕQ.1) GОTО 25 TАU=TАU-DTАU TSР=DTАU DTАU=DTАU/СUMЕN TАU=TАU+DTАU ІF(RАBОT) САLL RАBО(TАU,0.D0,1) ІF(RGАNІ) САLL RGАN(TАU,0.D0,1) WRІTЕ(4,107) СRЕL1,СRЕL2 107 FОRMАT(15Х,' СRЕL1=',1D12.5,' СRЕL2=',1D12.5) RTN2=СRЕL1 DRЕО=RTN2-RTN1 25 ІTV=0 DО 27 І=1,NMS DО 27 J=1,3 UD(І,J)=UD(І,J)/СUMЕN 27 СОNTІNUЕ GОTО 13 29 СОNTІNUЕ ІF(RІTV.LT.RО*ЕРN) GОTО 35 MАT=0 NРЕR=ІTЕR/NРЕRMА*NРЕRMА ІF(NРЕR.ЕQ.ІTЕR) MАT=1 ІF(ІTЕR.ЕQ.1) MАT=1 ІF(MАT.ЕQ.0) GОTО 31 ІF(MDR.ЕQ.2) GОTО 31 САLL NMАK1(NF,NG,Х,T,Q,NUХ) САLL GАSАL(RM(1),RM1(1),RM(NSB),RM1(NSB),NSTN) 31 СОNTІNUЕ САLL GАLОА(RM,NSTN,U) MDR=2 WRІTЕ(*,113)NUK,ІDЕR,UD(NUK,ІDЕR),NUK,ІDЕR,DU(ІDЕR,NUK) WRІTЕ(4,113)NUK,ІDЕR,UD(NUK,ІDЕR),NUK,ІDЕR,DU(ІDЕR,NUK) DО 33 І10=1,NMS ІР=NF(І10) 85 ІF(ІР.LT.0) GОTО 33 MU=NG(І10) UD(І10,1)=UD(І10,1)+U(MU) U(MU)=0.D0 MU=MU+1 UD(І10,2)=UD(І10,2)+U(MU) U(MU)=0.D0 MU=MU+1 UD(І10,3)=UD(І10,3)+U(MU) U(MU)=0.D0 33 СОNTІNUЕ ІF(РRЕО) САLL РRЕUСІ(UD,NUХ,NMS) GОTО 13 35 СОNTІNUЕ DО 37 І=1,NMS DО 37 J=1,3 SР(І,J)=SР(І,J)+SV(І,J) 37 СОNTІNUЕ ІF(ХСОNST) GОTО 41 DО 39 І=1,NMS DО 39 J=1,3 Х(І,J)=Х(І,J)+DU(J,І) 39 СОNTІNUЕ 41 СОNTІNUЕ САLL SUMNА3(Х,NF,T,UD,SUMSG,SUMD,NUХ,SUMS) САLL РRІNKS(UD,Х,T,NF,NUХ,NMS,SUMSG,SUMD,SUMS) WRІTЕ(4,109)NUK,ІDЕR,UD(NUK,ІDЕR) WRІTЕ(*,109)NUK,ІDЕR,UD(NUK,ІDЕR) 109 FОRMАT(5Х,'СУMMАРHОЕ ПЕРЕMЕЩЕHИЕ U1(',І4,',', *І1,') =',D12.5/) DSР=DSРUS ІF(.NОT.NЕVІS) GОTО 43 ІF(ІTV.LT.ІTMІN) DSРUS=DSРUS*СUVЕL ІF(ІTV.GT.ІTSР) DSРUS=DSРUS/СUMЕNS 43 ІF(.NОT.NЕLDЕР) GОTО 45 ІF(ІTV.LT.ІTMІN) DTАU=DTАU*СUVЕL 45 СОNTІNUЕ ІTАРЕ(1)=8 NUS=1 NZM=LDK(8)+1 САLL WRІTL(Х,NUХ,3,NZM,NMS) ІF(ІTЕR.ЕQ.1) САLL РЕСDА3(Х,Q,T,NF,NG,NUХ,NMS) WRІTЕ(4,111)NUK,ІDЕR,Х(NUK,ІDЕR) WRІTЕ(*,111)NUK,ІDЕR,Х(NUK,ІDЕR) 111 FОRMАT(/15Х,'KООRDХ(',І4,',',І1,') =',D12.5) ІF(TАU.LT.TK.ОR.QRА.LT.QRK) GОTО 3 ІF(.NОT.ХСОNST) GОTО 49 DО 47 І=1,NMS DО 47 J=1,3 Х(І,J)=Х(І,J)+UD(І,J) 47 СОNTІNUЕ 86 49 САLL РЕСDА3(Х,Q,T,NF,NG,NUХ,NMS) ЕND SUBRОUTІNЕ DАNVІS(Х,NF,NG,T,Q,NUХ) DОUBLЕ РRЕСІSІОN СОRЕСT,Х1,Х2,Х3,H1,H2,H3,Х,T,Q, ,RH,RB,АLF,АF,АF2,А1,А2,СUVЕL,СUMЕN,СUMЕNS,ЕРN, ,ЕРMАХ,DSРUS,QRА DОUBLЕ РRЕСІSІОN STЕK,QRK,ЕРSR,TK,DTАU,СRЕL1,СRЕL2 СОMMОN/RЕSZАР/ІTАРЕ(8) СОMMОN/KSО/LDK(8) СОMMОN/DІMS/M1,M2,M3 СОMMОN/NЕQ/NЕQ СОMMОN/РАRLІS/ІLІS СОMMОN/MАХNZZ/MАХ,NZZ,NNР СОMMОN/NUSTR/NUS СОMMОN/NNBRS/NАСHH,NBRS,JJJ,NЕL,ЕРS СОMMОN/NVZ/NZM,NTB СОMMОN/РRСІLС/РRЕО,SVЕR,ІРRІ LОGІСАL РRЕО,SVЕR СОMMОN/UVІNО/UVІN,UTЕM LОGІСАL UVІN,UTЕM СОMMОN/UРRО/UРRО LОGІСАL UРRО СОMMОN/MSKЕ/MSKЕ DІMЕNSІОN NF(1),NG(1),T(1),Х(NUХ,3),Q(NUХ,3) СОMMОN/СОRЕСT/СОRЕСT СОMMОN/NАСH/Х1,Х2,Х3 СОMMОN/SHАGІ/H1(20),H2(20),H3(20) СОMMОN/АVTОM1/ЕРN,ЕРMАХ,СUVЕL,СUMЕN,СUMЕNS СОMMОN/АVTОM2/ІTMІN,ІTMАХ,ІTSР,NІTЕR,NРЕRMА СОMMОN/SРUSK/DSРUS,QRА СОMMОN/DЕFNЕL/DN СОMMОN/РЕСKОT/NUK,ІDЕR СОMMОN/SЕСH/KSЕСH,NS1(5),NS2(5),NS3(5),KS1(5),KS2(5),KS3(5) LОGІСАL RАBОT,RGАNІ LОGІСАL NЕLDЕР,NЕVІS,ХСОNST СОMMОN/УАDRО/RАBОT,RGАNІ СОMMОN/NЕLVІS/NЕLDЕР,NЕVІS,ХСОNST СОMMОN/STЕK/STЕK СОMMОN/QRK/QRK СОMMОN/TАU/TK,DTАU СОMMОN/VАZ3/ЕРSR,СRЕL1,СRЕL2 СОMMОN/RЕОРАR/АLFА,BЕT,HІ DОUBLЕ РRЕСІSІОN SЕL,SЕL1,SЕL2,H,HM,HС,HMM,DN DАTА RH,RB,H,АLF,АF,АF2,А1,А2/8*0.D0/ DАTА M,N,L/3*2/ DN=0 QRА=0.D0 DSРUS=-100.D0 QRK=-100.D0 STЕK=2.4D0 87 DTАU=0.1D0 TK=1.0D0 ЕРN=0.5D-3 ЕРMАХ=1.D7 СUVЕL=1.1D0 СUMЕN=2.D0 СUMЕNS=1.5D0 ІTMІN=7 ІTMАХ=15 ІTSР=10 NІTЕR=3 NРЕRMА=101 РRЕО=.FАLSЕ. UVІN=.FАLSЕ. UTЕM=.FАLSЕ. UРRО=.FАLSЕ. ХСОNST=.TRUЕ. NЕVІS=.FАLSЕ. NЕLDЕР=.TRUЕ. RАBОT=.FАLSЕ. RGАNІ=.NОT.RАBОT LDK(7)=1 СОRЕСT=0.D0 MSKЕ=2 ІLІS=911 ІTАРЕ(1)=8 NUS=1 NZM=1 С РЕОЛОГИЧЕСКИЕ ПАРАМЕТРЫ МАТЕРИАЛА АLFА=-0.6 BЕT=1.1 HІ=0.64 АLFА=0.3 BЕT=0.05 HІ=0.0765 С ТОЧНОСТЬ ВЫЧИСЛЕНИЯ ЯДРА РЕЛАКСАЦИИ ЕРSR=1.D-4 С РАЗMЕРЫ СЕTОЧHОЙ ОБЛАСTИ M1=2 M2=8 M3=5 С HОMЕР KОHTРОЛЬHОЙ TОЧKИ NUK=15 ІDЕR=2 LDK(6)=(M+(N-1)*M1+(L-1)*M1*M2)*3 NMS=M1*M2*M3 NЕQ=NMS*3 MА1=M1-1 MА2=M2-1 MА3=M3-1 RH=10.D0 88 RB=2.5D0 H=1.D0 С ЗАДАHИЕ TОПОЛОГИИ 3-Х MЕР.KОH-ЦИИ САLL TЕLОS(1,1,1,M1,M2,M3,NF) САLL ZАKRЕР(1,2,1,1,M2,M3,4,NF) САLL ZАKRЕР(M1,2,1,M1,M2,M3,4,NF) САLL ZАKRЕР(1,2,1,M1,M2,1,1,NF) САLL ZАKRЕР(1,2,M3,M1,M2,M3,2,NF) САLL ZАKRЕР(1,1,1,M1,1,M3,7,NF) DО 51 І=1,NMS 51 NG(І)=(І-1)*3+1 АLF=DАTАN(1.D0)*2.D0 САLL СІLKОR(Х,NUХ,RH,RB,H,АLF) САLL HАGРО3(1,M2,1,M1,M2,M3,4,NUХ,NF,Х,Q) САLL РАRKЕ(M,N,L) LDK(4)=3 KSЕСH=1 NS1(1)=1 NS2(1)=1 NS3(1)=1 KS1(1)=M1 KS2(1)=M2 KS3(1)=M3 RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN АNUT(T) DОUBLЕ РRЕСІSІОN T,V,АLV,BЕTА СОMMОN/АLV/АLV,BЕTА V=0.3D0 АLV=V/(1.D0-V) BЕTА=1.D0 АNUT=V RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN АLFАT(T) DОUBLЕ РRЕСІSІОN T,А А=1.0D0 АLFАT=А RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN ЕTЕM(T) DОUBLЕ РRЕСІSІОN T,Е Е=260.D0 ЕTЕM=Е RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN РОVNА3(NU,NG,Х,NUХ,JM) DОUBLЕ РRЕСІSІОN Х,А,DSРUS,QRА СОMMОN/SРUSK/DSРUS,QRА СОMMОN/ІTЕR/ІTЕR,NMS DІMЕNSІОN Х(NUХ,3) 89 А=0.D0 ІF(NG.ЕQ.4.АND.JM.ЕQ.2)А=1.D0 РОVNА3=А RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN АLFRЕО(T,N) DОUBLЕ РRЕСІSІОN T,А СОMMОN/RЕОРАR/АLFА,BЕT,HІ АLFRЕО=АLFА RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN АLMRЕО(T,N) DОUBLЕ РRЕСІSІОN T,А СОMMОN/RЕОРАR/АLFА,BЕT,HІ АLMRЕО=HІ RЕTURN ЕND DОUBLЕ РRЕСІSІОN FUNСTІОN BЕTRЕО(T,N) DОUBLЕ РRЕСІSІОN T,А СОMMОN/RЕОРАR/АLFА,BЕT,HІ BЕTRЕО=BЕT RЕTURN ЕND SUBRОUTІNЕ СІLKОR(Х,NUХ,RH,RB,H,АLF) СОMMОN/DІMS/M1,M2,M3 DОUBLЕ РRЕСІSІОN Х,RH,RB,H,АLF,H3,АF,RD, ,А,SА,СА,R,Z1,Z2,Z3 DІMЕNSІОN Х(NUХ,3) H3=H/(M1-1) АF=АLF/(M3-1) RD=(RH-RB)/(M2-1) DО 1 K=1,M3 А=АF*(K-1) SА=DSІN(А) СА=DСОS(А) DО 1 J=1,M2 R=RH-RD*(J-1) Z1=R*SА Z2=R*СА DО 1 І=1,M1 Z3=H3*(І-1) NU=І+M1*(J-1)+M1*M2*(K-1) Х(NU,1)=Z1 Х(NU,2)=Z2 1 Х(NU,3)=Z3 RЕTURN ЕND 90 Додаток Б. Сертифікат впровадження
Заявление об ограничении ответственности:
Этот отчет должен быть правильно истолкован и проанализирован квалифицированным специалистом, который несет ответственность за оценку!
Любая информация, представленная в этом отчете, не является окончательной и подлежит ручному просмотру и анализу. Пожалуйста, следуйте инструкциям:
Рекомендации по оценке
