Проанализированный документ: Kozub_YuG.pdf Лицензия: ВОЛОДИМИР МАТІЄВСЬКИЙ
Детальный анализ тела документа:
Диаграмма соотношения частей:
Детали обработанных ресурсов:
250 - ОК /
12 - Ошибок
Активные ссылки (URL-адреса, извлеченные из документа):
Детальный анализ документа:
МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ ДЕРЖАВНИЙ ЗАКЛАД
«ЛУГАНСЬКИЙ НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ ІМЕНІ ТАРАСА ШЕВЧЕНКА»
Навчально-науковий інститут фізики, математики та інформаційних технологій Кафедра фізико-технічних систем та інформатики Козуб Юрій Гордійович МОДЕЛЮВАННЯ СКЛАДНИХ СИСТЕМ НА ОСНОВІ МЕТОДУ СКІНЧЕННИХ ЕЛЕМЕНТІВ. кваліфікаційна робота здобувача вищої освіти другого (магістерського) рівня освітньої програми
«Комп’ютерні науки та інформаційні технології»
за спеціальністю 122 Особистий підпис – ______________ Науковий керівник – _______ Смагіна О.О., к.пед.н., доцент (підпис) (посада, науковий ступінь, наукове звання, ініціали, прізвище) Полтава – 2023 АНОТАЦІЯ Козуб Ю.Г. Тема: Моделювання складних систем на основі методу скнченних. Спеціальність: 122 Установа: ДЗ ЛНУ імені Тараса Шевченка, 2023 р. Кваліфікаційна робота містить:80 стор., 29 рис., 48 джерел. Мета і задачі досліджень. Метою роботи є розробка ефективного методу моделювання просторових конструкцій із слабостисливих еластомерів та композитів з еластомерною матрицею Об’єкт дослідження – процес деформування елементів конструкцій з еластомерів та еластомерних композитів з урахуванням початкових напружень. Предмет дослідження – чисельне методи розв’язання варіацйних задач пружності конструкцій в умовах термосилового навантаження. Наукова новизна результатів роботи полягає у створенні на основі методу скінченних елементів нового ефективного підходу до чисельного моделювання зв’язаної задачі термов’язкопружності еластомерних та композитних конструкцій з початковими напруженнями під дією силових та температурних навантажень. При цьому: Ключові слова: метод скінченних елементів, НДС, термопружність, нестаціонарна теплопровідність, МІРЕЛА+ ЗМІСТ ВСТУП 4 РОЗДІЛ 1. ОГЛЯД СТАНУ ДОСЛІДЖЕНЬ З ПРОБЛЕМ МОДЕЛЮВАННЯ СКЛАДНИХ КОНСТРУКЦІЙ 6 1.1. Проблеми та задачі обчислювальної механіки складних тіл ...................... 6 1.2. Сучасні підходи до побудови моделей твердих тіл у тривимірній постановці ................................................................................................................ 8 РОЗДІЛ 2. МСЕ ДЛЯ ТРИВИМІРНИХ ЗАДАЧ ДЕФОРМУВАННЯ КОНСТРУКЦІЙ З ЕЛАСТОМЕРІВ ТА КОМПОЗИТІВ 17 2.1. Скінченний елемент для слабостисливих еластомерних конструкцій з початковими напруженнями ................................................................................ 17 2.2. Побудова матриці жорсткості для шаруватих конструкцій з урахуванням анізотропії властивостей ...................................................................................... 29 2.3. Побудова співвідношень скінченно-елементої моделі в’язкопружного деформування ........................................................................................................ 40 2.4. Дослідження збіжності розв’язків ................................................................ 45 РОЗДІЛ 3. ОСОБЛИВОСТІ ПРОГРАМНОЇ РЕАЛІЗАЦІЇ 52 РОЗДІЛ 4. НАПРУЖЕНО-ДЕФОРМОВАНИЙ СТАН ЕЛАСТОМЕРНИХ ТА ШАРУВАТИХ ЕЛЕМЕНТІВ АМОРТИЗАТОРІВ І ВІБРОСЕЙСМОІЗОЛЯТОРІВ 61 4.1 Дослідження амортизаторів складної форми для вібраційних машин ...... 62 4.2. Визначення параметрів віброізоляції конічних дробарок з гумометалевими елементами ............................................................................... 68 ВИСНОВКИ 75 СПИСОК ВИКОРИСТАНИХ ДЖЕРЕЛ 76 ВСТУП Актуальність теми. Розв’язання задач міцності, стійкості та довговічності для складних конструкцій можливе лише на основі чисельних методів, серед яких провідне місце займає метод скінченних елементів (МСЕ), що засвідчено його широким використанням в інженерній практиці. Проте питання дослідження на основі МСЕ конструкцій з еластомерів та еластомерних композитів не знайшли достатнього відображення в науковій літературі. Таким чином, розробка та реалізація об’єктно-орієнтованих ефективних підходів до розв’язання просторових зв’язаних задач термопружного деформування конструкцій з еластомерів і еластомерних композитів та дослідження на цій основі закономірностей їх термомеханічної поведінки є актуальною науковою проблемою сучасної обчислювальної механіки та представляє практичний інтерес. Мета і задачі досліджень. Метою роботи є розробка ефективного методу моделювання просторових конструкцій із слабостисливих еластомерів та композитів з еластомерною матрицеюю Об’єкт дослідження – процеси термопружного деформування елементів конструкцій з еластомерів та еластомерних композитів з урахуванням початкових напружень. Предмет дослідження – чисельне методи розв’язання варіацйних задач пружності конструкцій в умовах термосилового навантаження. Наукова новизна результатів роботи полягає у створенні
на основі методу скінченних елементів нового ефективного підходу до чисельного моделювання зв’язаної задачі термов’язкопружності еластомерних та композитних конструкцій з початковими напруженнями під дією силових та температурних навантажень. При цьому: Розроблено ефективну методику розв’язання задач
термов’язкопружності шаруватих конструкцій з анізотропними та ізотропними шарами. Побудовано матрицю жорсткості шаруватого суперелементу з урахуванням анізотропії шарів та початкових напружень в шарах конструкції. РОЗДІЛ 1. ОГЛЯД СТАНУ ДОСЛІДЖЕНЬ З ПРОБЛЕМ МОДЕЛЮВАННЯ СКЛАДНИХ КОНСТРУКЦІЙ 1.1. Проблеми та задачі обчислювальної механіки складних тіл Особливу увагу дослідників завжди привертали питання описання реакції матеріалів на механічне і теплове навантаження, вплив зовнішнього середовища на матеріал конструкції, а також аналіз умов і механізмів руйнування. З цією метою використовуються моделі поведінки матеріалів, які фактично представляють собою аналітичні вирази й рівняння, що характеризують реакцію матеріалу на відповідний вплив. В процесі експлуатації елементи конструкцій зазнають різного роду механічних і теплових навантажень, при яких поведінку матеріалу неможливо описати за допомогою тільки лінійних співвідношень. Більше того, вимоги зниження матеріаломісткості і собівартості конструкції, збільшення строків її роботи призводить до необхідності прогнозування поведінки матеріалу в нелінійній області. Обґрунтування основних рівнянь механіки взаємозв'язаних полів в матеріалах і елементах конструкцій, в тому числі і у неоднорідних висвітлено в роботах В.А. Баженова, [1] Григоренка А.Я. [24], О.І. Гуляра, О.С. Сахарова[2], В.І. Дирди [3], В.Г. Карнаухова [4], В.В. Киричевского [5], А.Д. Коваленка [6], Ю.М. Коляно [7] та ін. Широкий розвиток комп’ютерних систем став поштовхом до розвитку ефективних чисельних методів, таких як метод скінченних різниць (МСР), варіаційно-різницевий метод (ВРМ), метод прямих (МП), метод Монте-Карло (ММК), метод скінченних елементів (МСЕ), метод граничних елементів (МГЕ). В першу чергу з чисельних методів необхідно відзначити МСР, який отримав значне поширення при розв’язанні задач теплопровідності, так і для інших задач математичної фізики. В основу МСР покладено ідею про заміну диференціальних операторів деякими різницевими аналогами. Внаслідок приходять до розв’язання систем алгебраїчних рівнянь. Питання апроксимації, збіжність, точності і стійкості МСР розглянуто в багатьох роботах [8]. Варто відзначити, що до теперішнього часу розроблено низку універсальних програм, зорієнтованих на розв’язання задач теплопровідності для однорідних тіл, в основу яких покладено МСР. Метод прямих є поєднанням МСР за однією координатою і аналітичних розв’язків за рештою. Даний метод є ефективним у окремих практичних задачах для областей простої геометричної форми. В основі методу Монте-Карло лежить моделювання статистичного експерименту на комп‘ютері і реєстрація числових характеристик, отриманих з експерименту [9]. Серед варіаційних чисельних методів першим почали застосовувати варіаційно-різницевий метод. Головна ідея цього методу - заміна похідних під знаком інтегралу, який відповідає варіаційної постановці задачі. Цей метод з успіхом застосовувався для розв'язання складних задач теплопровідності, термопружності у роботах С. Г. Міхліна [9] та Ю. М. Шевченка [10]. Одним з найефективніших чисельних методів, що мають варіаційну основу, є метод скінченних елементів (МСЕ). Виникнення і розвиток МСЕ обумовлений у першу чергу роботам в галузі будівельної механіки. Застосування МСЕ дає можливість здійснювати розрахунки з урахуванням більш досконалих фізичних і математичних моделей. Для задач термопружності цей метод почав розвиватися порівняно нещодавно, на початку 60-х років минулого століття за кордоном і в 70-х роках на Україні. Великий вклад у його розвиток внесли вітчизняні та іноземні вчені В.А. Баженов [1], Д.В. Вайнберг [11], О.С. Городецький [12], О.І. Гуляр, О.С. Сахаров [21], В.В. Киричевський [5], Р. Галлагер [13], та ін. 1.2. Сучасні підходи до побудови моделей твердих тіл у тривимірній постановці При проектуванні конструкцій, що працюють в умовах високого теплового і механічного навантаження з високими техніко-економічними показниками, виникає необхідність прогнозування фізико-механічних характеристик і показників якості композиційних матеріалів у тривимірному аспекті. Дослідженню проблем і задач механіки неоднорідних ізотропних і анізотропних тіл присвячено значну кількість публікацій, серед яких особливе місце займають фундаментальні роботи А.Я. Александрова[14], О.М. Гузя [15], А.І. Лур’є [16], Б.Ю. Победрі [17], Ю.М. Работнова [18]. При вирішенні просторових і плоских задач широко використовується загальні подання рівнянь теорії пружності через гармонічні, бігармонічні, квазігармонічні і інші функції. А.І. Лур’є [16] надав загальний метод розв’язання задачі про рівновагу товстої плити, використовуючи запропонований ним символічний прийом розв’язання диференціальних рівнянь. Одним з універсальних і ефективних сучасних методів, як зазначалося раніше, при розв’язанні задач термопружності для анізотропних тіл зі складною геометрією є метод скінченних елементів. Для розв’язання просторових задач побудовано скінченні елементи різної форми - паралелепіпеди, тетраедри, трикутні призми та ін. Питанням комплексного розв’язання задач термопружності для однорідних та шаруватих тіл присвячено ряд праць. В роботі [19] проведено аналіз скінченно-елементних схем розрахунку багатошарових конструкцій і виведено наступні групи елементів: тривимірні, спеціалізовані двовимірні для розрахунку тришарових конструкцій, двовимірні для розрахунку шаруватих конструкцій з залежністю числа рівнянь від числа шарів, двовимірні для розрахунку шаруватих конструкцій з байдужістю розподільної системи до числа шарів. Останні в свою чергу поділяються на елементи, побудовані на основі гіпотези прямої нормалі, гіпотези прямою, більш загальних гіпотез, які враховують зсув і в деяких випадках обтиснення. У роботі Карнаухова В.Г. та ін. [20] розглянуто зв’язану тривимірну задачу про вимушені гармонічні коливання й дисипативний розігрів тришарової товстостінної в’язкопружної циліндричної панелі з шарнірним опиранням торців при дії на неї зовнішнього рівномірного поверхневого гармонічного тиску. Враховано взаємодію механічних і теплових полів за припущення про незалежність властивостей матеріалу від температури. Розроблено скінченно-елементний метод розв’язання динамічної задачі теплопровідності з відомим джерелом тепла. Розраховано амплітудно- й температурно-частотні характеристики, а також залежності власної частоти, коефіцієнта демпфірування, максимальних амплітуди й температури розігріву від товщини демпфуючого середнього шару панелі при її коливаннях на першій моді. При вирішенні зв'язаних задач термопружності для конструкцій зі анізотропних матеріалів застосовано різні теорії й підходи, які базуються на співвідношеннях зв'язаної задачі термов’язкопружності, отриманих О.А. Ільюшиним [21 116], Б.Є. Победрею [17], А.Д. Коваленком [6], В.Г. Карнауховим, І.К. Сенченковим, Я.О. Жуком [21] та ін. Одним з найважливіших критеріїв дослідження в'язкопружних тіл з нерівномірним температурним полем є врахування залежності від температури фізичних і теплофізичних характеристик матеріалу: модулю пружності або модулю зсуву, коефіцієнту температурного розширення й коефіцієнту теплопровідності. В цьому випадку вирішується зв'язана задача термов’язкоупружності. Розв'язання задач для багатошарових конструкцій складної форми пов'язано з математичними труднощами, які пов'язані з описом геометричних параметрів, умов теплообміну з зовнішнім середовищем, умов з'єднання, урахуванням внутрішніх джерел тепла, рішенням систем інтегральних рівнянь Для практичної реалізації є найбільш прийнятними уточнені теорії, що ґрунтується на безперервно-структурних моделях, в яких кількість і порядок дозвільних рівнянь не залежить від числа шарів. В роботі [22] розроблено скінченно-елементний метод розв’язування задачі про вимушені резонансні коливання й дисипативного розігріву в’язкопружних оболонок обертання з незалежним від температури властивостями з врахуванням поперечних деформацій. Як приклад, розраховано амплітудно- та температурно-частотні динамічні характеристики циліндричної панелі при коливаннях на першій згинній моді для випадків шарнірного та жорсткого опирання торців панелі. Досліджено вплив механічних граничних умов на вказані динамічні характеристики Для наближення розрахункової схеми конструкції до реальних умов необхідно враховувати перелічені чинники в їх взаємозв'язку. Вирішення такого складного завдання аналітичними методами пов'язано з непереборними труднощами. Для досягнення реальних успіхів при вирішенні зазначених задач доцільно використовувати чисельні методи, орієнтовані на використання сучасної обчислювальної техніки. Необхідно відзначити, що розрахунок шаруватих конструкцій можна здійснити на основі СЕ, побудованих для однорідних конструкцій. У цьому випадку кожен шар з різними фізико-механічними характеристиками піддається дискретизації на просторові (плоскі) елементи. Для двовимірних конструкцій з шарами порівнянної товщини, такий підхід використано Н.Н. Блумбергом [23] при дослідженні НДС багатошарових конструкцій при статичному і температурному навантаженні. Скінченно-елементна методика розв’язання просторової геометрично нелінійної задачі незв'язаної термопружності шаруватих анізотропних оболонок обертання при осесиметричному нагріванні і навантаженні застосовувалася Е.І. Гріголюком [25]. Проте при наявності в конструкції шарів малої товщини (клейовий прошарок, зовнішні шари тривимірної конструкції) застосування просторових СЕ складно, так як в цьому випадку потрібно забезпечити співмірність лінійних розмірів елементів. Для усунення такої невідповідності, розроблено різні типи СЕ, орієнтовані на розрахунок шаруватих композитних конструкцій, як специфічного класу. Проведений аналіз дозволяє виділити чотири основні групи СЕ, призначених для розрахунку даного класу конструкцій. До першої групи слід віднести СЕ, в яких існує залежність числа невідомих від кількості шарів, що призводить до розв’язання систем рівнянь високих порядків. Так, в роботі [26] запропоновано СЕ пластини, котрий у кожному вузлі має 2(n+1) ступенів свободи, що відповідають планарним (тангенціальних) переміщенням (n - число шарів) і одну ступінь свободи, відповідну нормальному переміщенню. Число ступенів свободи залежить від кількості шарів і незалежних переміщень в площині для жорстких шарів. Використовується варіант МСЕ - метод скінченних смуг. Число ступенів свободи для лінії вузлів з n жорсткими шарами одно 2(n+1). Криволінійний багатошаровий елемент для розрахунку вісесиметричних оболонок побудовано в [27]. Тут прийнято гіпотезу про сталість нормальних переміщень по товщині пакету, в жорстких шарах є незалежні планарні переміщення, а м'які шари чинять опір тільки поперечному зсуву, тобто кількість ступенів свободи СЕ залежить від кількості шарів. Другу групу складають спеціалізовані двовимірні СЕ побудовані для розрахунку тришарових композитних конструкцій. Так в роботі [28] розглядається два види трикутних СЕ для тришарової пластини. Для одного виду скінченних елементів побудова матриці жорсткості базується на певному законі розподілу компонент переміщення в окремих шарах елементів. Для іншого типу скінченних елементів при побудові матриці жорсткості базується на передбачуваному законі зміни напружень. Розрахунок тришарової пластини з легким заповнювачем на основі МСЕ розглянуто в роботі [29]. В основу покладено пропозицію про нехтування згинальної жорсткістю несучих шарів. Особливість полягає у низькому ступені поліномів, що апроксимують переміщення по області СЕ. Ряд інших типів СЕ для розрахунку тришарових композитних конструкцій розглянуто в роботі [30]. До третьої групи слід віднести двовимірні СЕ шаруватих композитних оболонок і пластин, побудовані на узагальненні гіпотезі Тимошенко-Рейсснера для всього пакету шарів, тобто гіпотези прямої. СЕ побудований для шаруватих конструкцій за аналогією з елементами однорідної конструкції. Наприклад, багатошарова анізотропна пластина з композитних матеріалів приводиться до квазіоднорідної шляхом використання гіпотез Тимошенко-Рейсснера для побудови ізопараметричних чотирикутних елементів. Анізотропні шаруваті системи на основі МСЕ розглянуто також у роботі [31]. Трикутний СЕ з дев'ятьома ступенями свободи на основі гіпотези Тимошенко-Рейсснера побудовано в [32]. Підхід, прийнятий при побудові третьої групи СЕ не дозволяє враховувати нелінійність зміни переміщення по товщині пакету. У зв'язку з викладеним, є перспективним підхід до побудови тривимірних СЕ на основі моделей, що враховують нелінійний закон зміни переміщень по товщині шарів. Ці елементи формують четверту групу. Так,в роботі [33] на основі моделі, що враховує деформації поперечного зсуву в шарах трансверсально-ізотропної плити, побудовано і досліджено прямокутний СЕ багатошарової ортотропної оболонки. У кожному вузлі призначено дев'ять ступенів свободи: по дві, зв'язані з планарними переміщеннями і функціями зсуву, і по три, що відповідають нормальним переміщенням і функціям обтиснення. Двовимірна схема розрахунку багатошарових трансверсально-ізотропних пластин, що враховує зміни нормальних переміщень від безпосереднього впливу зовнішніх навантажень розроблено В.Г. Пискуновим і представлено в роботі [34]. Значний вклад у розвиток підходів та методів механіки еластомерів і композиційних матеріалів та конструкцій з попереднім напруженням внесли вітчизняні вчені. Так можна відзначити роботи, виконані під керівництвом О.М. Гузя [35], Я.М. Григоренка [36]. Зазвичай нелінійне деформування конструкцій з переднапруженнями визначають двома шляхами. Перший шлях - експериментальний, другий - розрахунковий (аналітичний або чисельний розв’язок). Результати дослідження впливу початкового напруження на момент згину в пластині, що лежить на попередньо-напруженому напівпросторі, за умови дії рухомого навантаження наведено у статті С.Ю. Бабича, Ю.П. Глухова [37]. Рівняння руху пластини записано з врахуванням зсуву та інерції обертання. У просторі зображень одержано розв'язок в загальному випадку для стисливого і нестисливого напівпросторів, різних умов контакту та швидкостей руху навантаження. В рамках лінеаризованої теорії пружності для тіл з початковими напруженнями розглянуто постановку та метод розв’язку задач про поширення вісесиметричних хвиль в шаруватих композитних нестисливих заздалегідь напружених матеріалах при проковзуванні шарів. Досліджено випадок поширення хвиль вздовж шарів. Отримано дисперсійне рівняння для квазіпоперечних хвиль та його довгохвильові наближення. На основі лінеаризованої теорїї пружності отримано розв’язок вісесиметричної змішаної задачі про тиск пружного циліндричного штампу на шар з початковими (залишковими) напруженнями. Дослідження проведено в загальному вигляді для теорії великих початкових деформацій та різних варіантів теорії малих початкових деформацій при довільній структурі пружного потенціалу. У роботі [38] у формі гіперболічної системи диференціальних рівнянь першого порядку сформульовані визначальні рівняння руху пружного середовища із заданими початковими напруженнями. Виведено рівняння поширення малих збурень в попередньо напруженому ізотропному середовищі з довільною залежністю енергії пружною деформації від тензора деформацій і приведено рівняння для квадратичної залежності енергії пружної деформації від тензора деформацій. Одними з найбільш поширених композитних матеріалів є матеріали з в’язкопружною матрицею, серед яких важливе місце посідають еластомерні композити. Вивчення властивостей таких матеріалів, створення ефективних методів розрахунку та прогнозування ресурсу конструкцій з таких матеріалів є актуальним напрямком досліджень. Так у роботі [39] наведено результати натурних випробувань високонаповненого низькомодульного полімерного композиту з різними механічними складними видами впливу при різних температурах. Отримано нові дані про закономірності в'язкопружної поведінки композитних матеріалів у вигляді експериментальних залежностей в'язкопружних параметрів. В результаті дослідження закономірностей термомеханічної динамічної поведінки низько модульного високо наповненого полімеру отримано нові дані впливу умов навантаження на вязкоупружні параметри. Отримано нові дані про закономірності в'язкопружної поведінки композитних матеріалів у вигляді експериментальних залежностей в'язкопружних параметрів. В результаті дослідження закономірностей термомеханічної динамічної поведінки низькомодульного високонаповненого полімеру отримано нові дані впливу умов навантаження на в’язкопружні параметри. Розроблено феноменологічну багатофакторну математичну модель опису в'язкопружної поведінки конструкційного полімерного композитного матеріалу при складних видах впливу. В ході науково-дослідної роботи було розроблено основи експериментально-теоретичного апарату ідентифікації деформаційних властивостей в'язкопружних полімерів в умовах складних режимів навантаження з урахуванням можливостей новітніх випробувальних систем, засобів вимірювання і прикладного програмного забезпечення. Створено програму експериментальних досліджень динамічної
в'язкопружної поведінки матеріалу в умовах, близьких до експлуатаційних. Отримано верифіковані коефіцієнти багатофакторної феноменологічної математичної моделі опису в'язкопружної поведінки конструкційного полімерного композиту. Роботу [40] присвячено дослідженню пружних характеристик гумокордного матеріалу та створенню математичних моделей для дослідження напружено-деформованого стану конструкцій
із гумокордних матеріалів. Наведено підходи до визначення жорсткісних характеристик складних анізотропних гумокрдних матеріалів, які використовуються у каркасах та брекерах пневматичних шин. Проте при визначенні ефективних пружних характеристик гумокордного матеріалу поздовжний та поперечний модулі пружності корду вважалися рівними, що не для всіх матеріалів справедливо і що, відповідно, обмежує застосування наведених співвідношень. Математичну модель рівноважної конфігурації профілю пневматичної шини у вигляді інтегрального рівняння представлено в роботі [41]. У роботі [42] розглянуто задачу визначення деформацій повзучості лінійно-в’язкопружних матеріалів за різних режимів нестаціонарного навантаження, а саме ступеневого навантаження, повного розвантаження та циклічного навантаження. Розв’язок будується на основі спадкової теорії повзучості Больцмана – Вольтерра із дробово-експоненційним ядром. Нестаціонарні режими навантаження задаються за допомогою функцій Хевісайда. Результати розрахунків апробовано експериментально на задачах розрахунку деформацій нестаціонарної повзучості склопластиків, шаруватого пластику, полімербетону, міді, дуралюміну та нейлону. У роботі В. С. Саркісяна [43] на основі теорії нелінійної в'язкопружності побудовано розрахункову схему для визначення напружено-деформованого стану тонких пологих оболонок і пластин з урахуванням попереченого зсуву та нормального напруження. Розглянуті вище аналітичні і чисельні результати розв’язків задач термопружності, а також поділ типів СЕ багатошарових композитних конструкцій, показують, що останнім часом йде інтенсивна чисельна реалізація на основі МСЕ розроблюваних прикладних теорій. Однак кожна з методик має свою область застосування і не є універсальною. Найбільш ефективним є використання скінченно-елементної методики, побудованої на основі просторової постановки, оскільки в цьому випадку можна уникнути різних спрощуючих гіпотез і дозволити врахувати геометрію конструкції і просторовий характер термосилового навантаження. РОЗДІЛ 2. МСЕ ДЛЯ ТРИВИМІРНИХ ЗАДАЧ ДЕФОРМУВАННЯ КОНСТРУКЦІЙ З ЕЛАСТОМЕРІВ ТА КОМПОЗИТІВ 2.1. Скінченний елемент для слабостисливих еластомерних конструкцій з початковими напруженнями Формулювання вихідних положень почнемо з представлення процесу деформування у вигляді послідовності рівноважних станів ( ) ( ) ( ) ( ) ( ) n L 0 1 2 Ω ,Ω ,Ω ,.... Ω ,.... Ω , 0 L ( ) ( ) ( ) n де - начальний и кінцевий стан деформування відповідно; - , Ω Ω Ω довільний проміжний стан. ( ) ( ) ( ) 0 n n +1 Ω ,Ω ,Ω Рис. 2.1. Стани Вважається, що для всіх проміжних станів напруження, деформації, ( ) n переміщення відомі протягом всієї історії деформування до стану . Ω Завдання полягатиме в тому, щоб отримати рівняння для отримання ( ) ( ) n +1 n+1 змінних величин в стані в припущенні, що стан нескінченно Ω Ω ( ) n близько до стану і всі визначальні рівняння можна лінеаризувати по Ω відношенню до приросту змінних стану. Крок процесу деформування від стану ) ( ( ) n n +1 до стану можна (n+1) кроком. Ω Ω ( ) ( ) ( ) 1 0 + n n Ω Ω Ω , , Нехай положення довільної точки тіла в станах ( ) ( ) ( ) 0 n n+1 Q , Q , Q позначимо відповідно , радіуси-вектори цих точок 0 n n 1 + ( ) ( ) ( ) , , r r r . Тіло віднесено до базисної декартової системи координат і координати точок відповідно рівні (і = 1,2,3). ,Z , z Х і і і Радіус-вектори в такому випадку дорівнюють: 0 ( ) r = Х і , k k 0 n ( ) ( ) r = r u+ z= і Х = u і + ( )
k k k k k , n 1 0 + ( ) ( ) r = r u + u+ ∆Z =і z u= u +і + ∆ ( ) k k k k k k
, ( ) ( ) n n +1 де - вектори переміщень в станах , , відповідно. u, u + ∆u Ω Ω 0 0 ( ) ( ) n+1 n ε ε + ε Позначимо тензори деформацій Гріна в станах та через і Ω Ω іj іj іj відповідно: 1 1 n n 0 0 ( ) ( ) ( ) ( ) 0 ε =
r r r −r u= u u +u + ( ) іj і , j j ,і k ,і k , j ,і , j ,і , j ( ) 2 2 1 n 1 n 1 0 0 + + ( ) ( ) ( ) ( ) 0 ε + ε = r r − r r
= іj іj і j і j , , , , ( ) 2 1 ) ( ) ( ) ( ) ( ) ( = u + ∆u + u + ∆u + u + ∆u u + ∆u . (2.1) і і j j k k k k , j ,і , j ,і 2 ( ) n З іншого боку приймаючи за початковий стан і використовуючи Ω прямокутні координати приросту деформацій можна визначити n 1 n 1 n n + + ( ) ( ) ( ) ( ) u ∂∆ u u u ∂∆ ∂∆ ∂∆ 1 1 ∂r ∂r ∂r ∂r j * і k k . (2.2) ε = − = + + іj 2 z z z z 2 z z z z ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ і j і j j і і j Прирости деформацій пов’язані співвідношеннями: ∂Х ∂Х ∂z ∂z * * m n m n ε = ε ε = ε (2.3) іj mn , іj mn ∂z ∂z ∂Х ∂Х і j і j Компоненти тензора приросту деформації можна представити у вигляді лінійної і нелінійної складових: Л Н u ∂∆ u u u ∂∆ ∂∆ ∂∆ 1 1 j * * * і k k + . ε = ε + ε = + іj іj іj z z ∂z ∂z 2 ∂ ∂ 2 j j і і ( ) ( ) n n +1 . Якобіани перетворення = D z , D Z = і , j j , і ( ) = δ + + ∆ Z u u . = δ + z u , і , j іj і і і , j іj і , j , j Для описання напруженого стану введемо тензор напружень Ейлера. У іj ( ) n Q компоненти тензора напружень мають значення , що діють на σ точках 0 іj іj z − соnst ( z + dz ) − соnst гранях , , що діють на гранях , та σ + σ Z − соnst і і і і 0 нескінченно малого паралелепіпеду, що включає точку ( ) Z + dZ − соnst і і ( ) n +1 Q . ( ) n +1 Q У точці напруження можна описати модифікованим тензором Кирхгофа [44]. ∂z ∂z 1 j іj іj* kl kl і ( ) (2.4) σ + σ = σ + σ , 0 0 ( ) n ∂ ∂ Х Х D k l ∂z ∂z 1 j іj* kl і σ = σ . (2.5) n ( ) ∂ ∂ Х Х k l D Bаріаційна постановка задачі для пружного тіла з початковими напруженнями зі заданими додатковими масовими та поверхневими силами у ( ) n +1 формі принципу віртуальної роботи у стані має вигляд: Ω * n іj ( ) іj * 0 0 , (2.6) σ + σ δε + + δ − + δ∆ = q q u dV р р u dS 0 ( ) ( ) іj і і і і і і ( ) 0 ∫∫∫ ∫∫ V S де р , q – вектори поверхневих та об’ємних додаткових сил на n+1 кроці. і і ( ) n – врівноважений стан, то у рівнянні члени варіації пружної Якщо Ω енергії та з урахуванням початкового напруження іj 0 0 * n n [ ] q u dV р u dS 0 σ δε − δ − δ = . (2.7) ∫∫∫ ∫∫ іj і і і і 0 V S рівняння рівноваги приймає вигляд: u u ∂∆ ∂∆ 1 іj* іj ( ) ( ) n n * k k (2.8) dS . р u u dV q − δ∆ = − δ∆ σ δε − σ δ 0 ∫∫∫ ∫∫ і і іj і і 0 z z ∂ 2 ∂ j і V S Для дослідження напружено-деформованого стану просторових конструкцій з еластомерів розглянемо ізопараметричний скінченний елемент (СЕ) серендіпова сімейства у вигляді шестигранного паралелепіпеда з довжиною ребер, що дорівнює 2. Початок базисної системи координат Z і і довільної місцевої системи х , осі якої збігаються з напрямом його ребер, і поміщаємо в центрі куба. Розглянемо побудову матриць жорсткості скінченного елемента з початковими напруженнями на основі моментної схеми скінченних елементів. Під попередніми напруженнями розуміємо ті напруження, які виникли в конструкції в початковому стані, тобто перед початком процесу деформації, що розглядається, до прикладення робочого навантаження. х1 24 4 32 12 х3 0 1 21 х2 9 29 Рис. 2.2. Кубічний скінченний елемент У задачі з попередніми напруженнями виберемо початковий стан як відліковий. Переміщення за об’ємом СЕ серендіпова сімейства апроксимуємо у вигляді: lmn ( ) lmn ( ) рqr ψ = ω u (2.9) ∑ k k рqr ( ) рqr де - коефіцієнти розкладання, - набір степеневих функцій виду: ψ q р r х х х ( ) рqr 3 1 2 , = ψ (2.10) r! q! р! де (р=0,1,2,…,M; q=0,1,2,…,N; r=0,1,2,…,L) - степені апроксимації за відповідними напрямами координат. 000 100 010 001 110 011 ( ) ( ) ( ) ( ) ( ) ( ) = ω + ω + ω + ω + ω + ω + u
8%
8%
7%
6%
5%
5%
4%
4%
4%
3%
3%
3%
2%
2%
2%
2%
2%
1%
0,9%
0,9%
0,7%
0,7%
0,6%
0,6%
0,5%
0,5%
0,5%
0,4%
0,4%
0,3%
0,3%
0,3%
0,3%
0,3%
0,2%
0,1%
х х х х х х х 1 2 3 1 2 2 3 k k k k k k k 1 1 2 2 101 111 200 210 ( ) ( ) ( ) ( ) + ω + +ω + ω + ω + х х х х х х х х 1 1 3 2 3 1 1 2 k k k k 2 2 1 1 1 1 2 2 2 2 201 211 020 120 ( ) ( ) ( ) ( ) + ω + ω + ω + ω + х х х х х х х х 1 3 1 2 3 2 1 2 k k k k 2 2 2 2 1 1 1 1 2 2 2 2 021 121 002 102 ( ) ( ) ( ) ( ) + ω + ω + ω + ω + х х х х х х х х 2 3 1 2 3 3 1 3 k k k k 2 2 2 2 1 1 1 1 2 2 3 3 012 112 300 310 ( ) ( ) ( ) ( ) + ω + ω + ω + ω + х х х х х х х х 2 3 1 2 3 1 1 2 k k k k 2 2 6 6 1 1 1 1 3 3 3 3 301 311 030 130 ( ) ( ) ( ) ( ) + ω + ω + ω + ω + х х х х х х х х 1 3 1 2 3 2 1 2 k k k k
6 6 6 6
1 1 1 1 3 3 3 3 031 131 003 103 ( ) ( ) ( ) ( ) + ω + ω + ω + ω + х х х х х х х х 2 3 1 2 3 3 1 3 k k k k 6 6 6 6 1 1 3 3 013 113 ( ) ( ) + ω + ω х х х х х 2 3 1 2 3
k k 6 6 ( ) рqr Для координатних функцій справедливе співвідношення ψ диференціювання: ( ) ( ) ( ) α+β+γ рqr р −α q −β r −γ ∂ ψ = ψ , (2.11) β α γ де , , - порядок похідних. На основі розподілу переміщень отримуємо розподіл похідних функції переміщень: (100) (200) (110) (101) (111) (210) = ω + ω + ω + ω + ω + ω + u
8%
8%
7%
6%
5%
5%
4%
4%
4%
3%
3%
3%
2%
2%
2%
1%
0,9%
0,9%
0,7%
0,6%
0,5%
0,5%
0,5%
0,4%
0,3%
0,3%
0,2%
х х х х х х х k ,1 1 2 3 2 3 1 2 k k k k k k 1 1 (201) (211) (120) (102) 1 3 1 2 3 2 2 +ω + ω + ω + ω + х х х х х х х 2 3 k k k k 2 2 1 1 1 1 (121) (112) (300) (310) 2 2 2 2 + ω + ω + ω + ω + х х х х х х х 2 3 2 3 1 1 2 k k k k 2 2 2 2 1 1 1 1 (301) (311) (130) (131) 2 2 3 3 + ω + ω х х х х х + ω х + ω х х + 1 3 1 2 3 2 2 3 k k k k
6 6 2 6 1 1 (103) (113) 3 3 + ω х + ω х х , 3 2 3 k k 6 6 (010) (110) (020) (011) (111) (120) = ω + ω + ω + ω + ω + ω + u
8%
8%
7%
6%
5%
5%
4%
4%
4%
3%
3%
3%
2%
2%
2%
1%
0,9%
0,9%
0,7%
0,7%
0,5%
0,5%
0,5%
0,4%
х х х х х х х k ,2 1 2 3 1 3 1 2 k k k k k k 1 1 (121) (210) (021) (012) 2 2 +ω + ω + ω + ω + х х х х х х х 1 1 2 3 2 3 3 k k k k 2 2 1 1 1 1 (211) (112) (310) (013) 2 2 3 3 + ω + ω + ω + ω + х х х х х х х 1 3 1 3 1 3 2 k k k k 2 2 2 2 1 1 1 1 (031) (130) (131) (130) 2 3 2 2 + ω + ω х х х х + ω х + ω х х х + 2 3 2 1 2 3 1 32 k k k k
6 6 2 6 1 1 (311) (113) 3 3 + ω х х + ω х х ,
8%
8%
7%
6%
5%
5%
4%
4%
4%
3%
3%
3%
2%
2%
2%
2%
1%
0,9%
0,9%
0,7%
0,7%
0,6%
0,5%
0,5%
0,5%
0,4%
0,4%
0,3%
0,3%
0,2%
0,2%
0,1%
1 3 1 3 k k 6 6 (001) (101) (011) (002) (111) (012) = ω + ω + ω + ω + ω + ω + u х х х х х х х k ,3 1 2 3 1 2 2 3 k k k k k k 1 1 (102) (112) (201) (211) 1 3 1 2 3 2 2 +ω + ω + ω + ω + х х х х х х х х 1 1 2 k k k k 2 2 1 1 1 1 (021) (121) (103) (013) 2 2 2 2 + ω + ω + ω + ω + х х х х х х х 2 1 2 1 3 2 3 k k k k 2 2 2 2 1 1 1 1 (301) (311) (031) (113) 3 3 3 2 + ω + х х х ω х + ω х х + ω х + 1 1 2 2 1 2 3 k k k k 6 6 6 2 1 1 1
(131) (003) 3 2 + ω х х + ω х . 1 2 3 k k 6 2 енергії пружної деформації для СЕ У відповідності з (2.8) варіацію W δ приймаємо у вигляді: 1 іj іj ( ) (2.12) u dv. W u δ = σ δε + σ δ ∫∫∫ , j ,і k іj k 0 2 V Компоненти тензора пружних напружень приймаємо у вигляді закона Гука: іj іk іl іj σ = 2µg g ε + λθg , (2.13) kl µ ,λ - коефіцієнти Ляме; де іj g - компоненти метричного тензора; - функція зміни об’єму. θ Представимо компоненти тензора деформацій в вигляді розкладання ε іj в ряд Маклорена: L M N іj іj іj stg stg ( ) ( ) . е ε = ψ (2.14) ∑ ∑ ∑ іj іj 0 0 0 S t g = = = Вираз (2.13) в матричній формі буде мати вигляд: T { } { } { } ε = е ψ . (2.15) іj іj іj тензора деформації в ряд Маклорена Компоненти розкладання е ε іj іj ( ) рqr за допомогою перетворення: представимо через коефіцієнти ω k k е = F ω . (2.16) { } { } іj іj k Для урахування слабкої стисливості еластомерів функцію змінення об’єму представимо також у вигляді ряду [134]: 1 1 1 − − − l m n ( ) ( ) αβγ αβγ θ = ξ ψ , (2.17) ∑ ∑ ∑ 0 α = β= 0 0 γ= ( ) αβγ визначаються співвідношенням: де коефіцієнти розкладання ξ ( ) α+β+λ іj ( ) αβγ ∂ u g і , j (2.18) ξ = , β γ α ∂ ∂ ∂ х х х 1 2 3 х = х = х = 0 1 2 3 які у матричній формі мають вигляд: T { } {} (2.19) θ = ξ ψ ; ( ) θ k [ ] {} { }. ξ = F ω (2.20) θ k Підставляючи (2.13), (2.15) у варіацію енергії пружної деформації (2.12): T T іk jl { } } { } { { } 2 W е g g е dV δ = δ ψ µ ψ + ∫∫∫ ( ) ) (іj kl kl V T T іj + { } { } { } { } δξ g dV ψ λ ψ ξ + ∫∫∫ ( ) ( ) θ θ V 1 іj [ ] + σ δ + δ u u u u dV , ∫∫∫ k ,і k , j k ,і k , j 0 2 V T T T T іj θ іjkl W е Е е Е , δ = δ + δξ ξ + δω ψ σ ψ ω { } { } { } { } { } { } { } { } { } , , іj kl k k l k j k 0 1 T іjkl іk іl [ ] } { } { = µ ψ ψ 2 Е g g gdх dх dх ; де ∫∫ ∫ ( ) ( ) 1 2 3 іj kl −1 1 T іj θ [ ] { } { } Е g gdх dх dх ; = λ ψ ψ (2.21) ∫∫∫ ( ) ( ) 1 2 3 θ θ −1 З урахуванням (2.14) и (2.18) вираз буде мати вигляд: δW T T T T θ ( ) s іjkl t s t W F Е F F Е F δ =δ ω ω +δ ω ω + { } { } { } { } s іj kl t s t θ θ ( ) ( ) T T іj +δ ω ψ σ ψ ω . (2.22) { } { } { } { } { } k k , , k і k j 0 ( ) ( ) Для побудови матриці жорсткості необхідно перейти у виразі (2.22) від коефіцієнтів до значень переміщень у вузлах скінченного елементу, що ω s [ ] функцій форми можна здійснити за допомогою матриці перетворення . А N L Ця матриця задає зв'язок між апроксимуючими функціями Лагранжа і степеневими ψ: T s u = u N , (2.23) { } { } і і s ( ) де – функції форми; N ( ) s s u – вузлові значення переміщень СЕ. і З іншого боку: T { } { } u = ω ψ ′ ′ k k Між степеневими функціями та функціями форми існує зв’язок: T [ ] { } { } N А , = ψ (2.24) [ ]{ } { } ω = А u , k k де Аматриця переходу від степеневих функцій до функції форми. З урахуванням (2.22) вираз (2.21) набуде вигляду: T T T s іjkl t [ ] [ ][ ] [ ] [ ] { } { } δW =δ u А F Е F А u + s іj kl t T T T s t [ ] [ ] [ ] [ ] [ ] { } { } u А F Е F А u +δ + s θ θ θ t іj T T T [ ] [ ] [ ] [ ] [ ] { } { } + δ ψ σ ψ = u А А u t , j t s s ,і 0 T T st st [ ] [ ] { } { } { } { } =δ u K u + δ u K u (2.25) s t s t 0 T T T T st s іjkl t s t K А F Е F А А F Е F А = + де (2.26) [ ] [ ] [ ] [ ] [ ] іj kl θ θ θ ( ) ( ) 96 96 × – матриця жорсткості. Матриця впливу початкових напружень обчислюється за формулою: T T іj st s t . (2.27) = K А F H F А ] ] [ [ 0 (, ) (, ) k і k j 0 . Таким чином, система лінійних алгебраїчних рівнянь буде мати вигляд: s t s t s { } [ ] { } K + K u = Р . 0 t З виразів (2.25) и (2.26) видно, що для побудування матриць жорсткості st st s s [ ] K та необхідно визначити спеціальні матриці та . [ ] А F K F , k j 0 (, ) іj ( ) рqr Коефіцієнти розкладання (2.15) визначимо згідно МССЕ за е іj формулами [5]: рqr ( ) рqr ( ) 1v s µ+ η е b ; = ω ∑ s 11 ) ( 1 р q r + −µ −ν −η µνη рqr ( ) рqr ) (µν 1η s + е = ω b ; ∑ s 22 ( ) −η р −µq +1−νr µνη рqr ( ) рqr ( ) µνη+1 s е = ω b ; ∑ s 33 ( ) р −µq −νr +1−η µνη рqr 1 ( ) рqr ( ) ( ) µ+ νη + η µν s s 1 1 ) ( = ω + ω b ; е b ∑ s s 12 ( ) ( ) −µ −ν+ −η −µ+ −ν −η р q r р q r 1 1 2 µνη рqr 1 ( ) рqr ( ) ( ) s s 1 1 µνη+ µ+ νη ) ( е b b ; = ω + ω ∑ s s 13 ( ) ( ) р q r р q r 1 1 −µ+ −ν −η −µ −ν −η+ 2 µνη рqr 1 ( ) рqr ( ) ( ) µνη+ µν+ η 1 1 s s (2.28) ) ( = ω + ω е b b , ∑ s s 23 ( ) ( ) −µ −ν+ −η −µ −ν −η+ 1 1 р q r р q r 2 µνη ( ) µ+ν+η ∂ z s s = b . де ( ) µνη µ µ ν ∂ ∂ ∂ х х х 1 2 3 = = = х х х 0 1 2 3 ( ) рqr ( ) рqr е Аналіз кожної з компонент показав, що деякі коефіцієнти не ω іj k входять до розкладання для апроксимації переміщень (2.10). Ті коефіцієнти ( ) рqr е деформацій , які містять хоча б один з членів, відсутній в (2.10), повинні іj бути опущені в розкладанні (2.14). Для отримуємо наступне розкладання: ε 11 1 ( ) ( ) ( ) ( ) ( ) 000 010 110 110 020 2 ε ≈ е + е х + е х + е х х + е х + 11
1 2 1 2 11 11 11 11 11 2 2 ( ) ( ) ( ) ( ) 001 101 011 111 (2.29) + + + + + е х е х х е х х е х х х 3 1 3 2 3 1 2 3 11 11 11 11
1 1 1 ( ) ( ) ( ) 021 002 012 2 3 2 2 + + + х х х е х е х е 2 11 2 11 3 3 11 2 2 2 k F . На основі виразу (2.28) обраховується матриця 11 k k k k k F F F F F Аналогічно отримаємо матриці , , , , , 22 33 21 31 32 k k k F F F , , відповідно для тензорів та ε ε ε ε ε , , , , s
s s (,1) (,2) (,3) 22 33 21 31 12 u u u , , векторів похідних s s s ,1 ,2
,3 А Розглянемо процес побудови матриці перетворення . [ ] ~ Апроксимацію компонентів переміщень за об’ємом СЕ можна u і ( ) рqr u представити через функцію форми і вузлові переміщення : N ( ) рqr і lmn ( ) рqr ~ u = u N , ∑ ( ) і рqr і рqr У матричній формі: T ~ (2.30) { } { } = u u N , і і T } } { { u = ω ψ k k Як видно з (2.30), за допомогою лінійного перетворення можна перейти від системи функцій форми до системи координатних функцій шляхом співвідношення: T . (2.31) { } { } [ ] ψ N = А Матриця визначається з співвідношень (2.30). [ ] А 2.2. Побудова матриці жорсткості для шаруватих конструкцій з урахуванням анізотропії властивостей Конструкції з композиційних матеріалів знаходять широке застосування у різних галузях народного господарства. Використання конструкцій з полімерних композиційних матеріалів дозволяє знизити матеріаломісткість конструкцій, скоротити терміни виробництва, підвищити корозійну стійкість виробів. Одними з широко використовуваних матеріалів є шаруваті композиційні матеріали з еластомерною матрицею, що складаються з тонких односпрямованих армованих шарів або з почергових односпрямовано армованих та ізотропних шарів. До того ж кількість таких шарів в пакеті матеріалу може сягати двозначних величин. Зокрема, великогабаритні шини мають до 40 шарів гуми та гумокорду на боковині та ще більше на біговій частині, якщо враховувати шари брекеру. Застосування шаруватих елементів конструкцій обумовлене їх більшою ефективністю в порівнянні з однорідними деталями. За допомогою шаруватості досягається оптимальне співвідношення жорсткісних і демпфуючих характеристик багатьох елементів віброізолюючих систем. У сучасних промислових конструкціях, виготовлених з шаруватих композитних матеріалів, при експлуатаційних навантаженнях відбувається значний нагрів. При цьому істотний вплив на характер напружено- деформованого стану дає залежність фізико-механічних характеристик від температури, що приводить до суттєвих змін форми та в ряді випадків до руйнування. В загальному випадку необхідно вирішувати зв’язану задачу термопружності. Для скінченно-елементного розв’язання задачі термопружності для просторових шаруватих конструкцій за умови дії попередніх навантажень, що мають складну геометричну форму можна використовувати різні підходи. Розглянемо суперелементний підхід до моделювання поведінки шаруватої конструкції. Рис. 2.3. Шаруватий скінченний елемент У загальному випадку такий елемент може бути складений з n шарів (елементів). Кожен з таких шарів має або анізотропні властивості, або ізотропні властивості слабкостисливого еластомеру. Опис напружено-деформованого стану композиційних шарів здійснюється за допомогою співвідношень, що враховують специфічні особливості анізотропних матеріалів. Ці особливості матеріалу можна задавати фізичними співвідношеннями, що зв'язують напруження і деформації, у вигляді закону Гука: іj іjkl σ = ε С , (2.32) kl іj - компоненти тензора напружень; де σ іjkl – компоненти тензора пружних констант матеріалу; С – компоненти тензору деформацій. ε kl Тензор пружних сталих анізотропного матеріалу має властивість симетрії: іjkl jіkl іjlk іjkl іjkl lkіj . (2.33) С = С = С С ; С ; С = ∗ ∗ ∗ ∗ ∗ ∗ та Для опису скінченного елементу використовується базисна зв’язана зі скінченним елементом місцева ортогональні системи координат. Місцева система координат починає відлік в центрі k –го шару скінченного елемента. Між базисними та місцевими координатами елементу, що моделює поведінку шару, прийнято відображення: 8 1 2 3 і l і (2.34) ( ) = z N х , х , х z , ∑ l 1 = l 1 l 1 1 2 2 3 3 і ( )( )( ) N = 1 + х х 1 + х х 1 + х х z де - функції форми для -го вузла СЕ, - l l l l 8 глобальні координати l-го вузла СЕ. Для опису переміщень шаруватого композиту приймемо гіпотезу про нерозривність поля переміщень на межі поділу шарів: ( ) ( ) k + k +1 − . (2.35) u = u Компоненти вектору переміщень точок k-го шару СЕ можна також апроксимувати у вигляді розкладання ) ( ) ( ) ( ) (000 100 010 001 ( ) ( ) ( ) 100 010 001 = ω + ω ψ + ω ψ + ω ψ + u
і і і і і (2.36) ( ) ( ) ( ) ( ) 110 101 011 111 ( ) ( ) ( ) ( ) 110 101 011 111 + ω ψ + ω ψ + ω ψ + ω ψ , і і і і
1 1 1 ( ) рqr ( ) рqr u = ω ψ , або (2.37) ∑ ∑ ∑ і і 0 0 0 р = q = r = ( ) рqr ω де - коефіцієнти розкладання; і ( ) рqr - набір степеневих координатних функцій виду: ψ р q r х х х ( ) рqr 3 1 2 ( ) , р ,q , r , . 0 1 ψ = = (2.38) р! q! r! ( ) рqr виконується умова диференціювання: Для координатних функцій ψ ( ) ( ) ( ) α +β+ γ ρ− α −β − γ рqr q r . (2.39) ∂ ψ = ψ Похідні функції переміщень в середині k-го шару мають вигляд: ) ( ) ( ) ( ) (001 111 100 110 ( ( ) ( ) ) 010 001 011 , u = ω + ω ψ + ω ψ + ω ψ 1
і , і і і і ( ) ( ) ( ) ( ) 010 110 011 111 ( ) ( ) ( ) 100 001 101 , u = ω + ω ψ + ω ψ + ω ψ і ,2 і і і і ( ) ( ) ( ) ( ) 001 101 011 111 ( ) ( ) ( ) 100 010 110 . (2.40) u = ω + ω ψ + ω ψ + ω ψ 3 і , і і і і
) ( рqr та Між компонентами розкладання переміщення по базисних функціях ψ ( ) l існує взаємний зв’язок: переміщеннями вузлових точок шару u і T ) (l T { } } { { } { } , (2.41) u u N = = ω ψ і і і k ( ) ω = А u де ; (2.42) { } [ ] і і { } T . (2.43) [ ] { } { } = ψ А N Компоненти тензора деформацій також представимо у вигляді розкладання по степеневим функціям: ( ) stq ( ) stq . (2.44) ε = ψ е ∑ іj іj stq В матричній формі: T , (2.45) { } { } е ψ = ε ( ) іj іj іj ∂ ∂ u u 1 j і k k де - лінійний тензор деформації; ε = + С С іj k k j і 2 ∂х ∂х і ∂z і С – тензор перетворення координат. k k ∂х Представимо коефіцієнти розкладання тензора деформацій , що е ε іj іj ) ( рqr ω (2.36): входять в (2.45) через коефіцієнти апроксимації переміщень і s . (2.46) е = F ω { } { } іj іj s S ] [ можна отримати у відповідності до моментної схеми Матрицю F іj скінченних елементів (2.27) [233]: ( ) рqr Аналізуючи кожну з компонент помічаємо, що деякі коефіцієнти е 11 ( ) mnl ω не входять в розкладання для апроксимації переміщень (2.43). s ( ) рqr Коефіцієнти розкладання деформацій , які містять хоча б один з членів, е 11 що відсутні в (2.40), мають бути вилучені в розкладанні (2.44). Таким чином отримуємо наступне розкладання: ( ) ( ) ( ) ( ) 000 001 010 011 ( ) ( ) ( ) 001 010 011 . (2.47) е е е е ε = + ψ + ψ + ψ 11 11 11 11 11 s F . На основі розкладання (2.46), (2.47) будується матриця 11 Проводячи аналогічно викладки для інших компонент тензора деформацій, отримуємо наступні вирази їх розкладання в ряд: ) ( ( ) ) ) ( (000 100 001 101 ( ) ( ) ) (100 001 101 , ε = е + е ψ + е ψ + е ψ 22 22 22 22 22 ( ) ( ) ( ) ( ) 000 100 010 110 ( ) ( ) ( ) 100 010 110 , ε = е + е ψ + е ψ + е ψ 33 33 33 33 33 ) ( ) ( 000 001 ) (001 , е е ε = + ψ 12 12 12 ( ) ) (000 010 ) (010 , е е ε = + ψ 13 13 13 ( ) ( ) 000 100 ( ) 100 . (2.48) ε = е + е ψ 23 23 23 s s s s s F F F F F , , , , будуються аналогічно. Матриці 22 33 23 12 13 Варіація пружної енергії деформації анізотропного шару з початковими напруженнями має вигляд: 1 іj іj ( ) W dV u u dV , δ = σ δε + σ δ ∫∫∫ ∫∫∫ іj k ,і r , j 0 2 V V 1 іj T іjkl { } [ ] { } С δW = δε ε dV u u u u dV . + σ δ + δ (2.49) ∫∫∫ ∫∫∫ іj kl k ,і k , j k ,і k , j 0 2 V V Підставляючи вираз (2.44) в (2.49) отримуємо T T T іj іjkl δW = е δ ψ С ψ е dv + δω ψ σ ψ ω = { } { } { } { } { } { } { } { } { } ( , ) ( , ) іj kl k k і k j k іj kl 0 ∫∫∫ ( ) ( ) V T T T s іjkl t іj = δω F H F ω + δω H ω = { } { } { } { } s іj kl t k k T s іjkl t = δu А F H F А u + { } } { [ ] [ ] s іj kl t T T T іj s t , (2.50) +δ u А F H F А u { } { } [ ] [ ] s k (,і ) k (, j ) t 0 1 1 1 T іjnm іjnm [ ] [ ]{ } { } = ψ ψ dх dх dх H С g тут ∫ ∫ ∫ ( ) ( ) 3 2 1 іj mn −
1−1−1 1 1 1 T іj іj = ψ σ ψ H gdх dх dх (2.51) { } { } 0 k і ( , ) 1 2 3 k j , 0 ∫ ∫ ∫ ( ) − − − 1 1 1
- матриця пружних констант, що характеризують фізичні властивості анізотропного матеріалу k-го шару та матриця впливу початкових напружень. іjn m [ ] H Для побудови матриці компоненти тензора пружних постійних іjkl задаються відповідно до виду анізотропії композиту в системі координат С * пов'язаної з осями анізотропії. ξ , і Матриця пружних сталих ортотропного тіла в цій системі координат має вигляд: 1111 1122 1133 С С С
8%
6%
3%
2%
2%
2%
0,7%
0,7%
0,2%
0 0 0 0 0 0 * * * 1212 1221 С С 0 0 0 0 0 0 0 * * 1313 1331 С С 0 0 0 0 0 0 0 * * 2112 2121 С С 0 0 0 0 0 0 0 * * іjkl 2211 2222 2233 [ ] С = С С С 0 0 0 0 0 0 * * * * 2323 2332 С С 0 0 0 0 0 0 0 * * 3113 3131 С С 0 0 0 0 0 0 0 * * 3223 3232 С С 0 0 0 0 0 0 0 * * 3311 3322 3333 С С С 0 0 0 0 0 0
* * * Компоненти тензора пружних сталих ортотропного матеріалу обчислюються за допомогою співвідношень: Е Е 1122 1111 2 1 ( ) ( ); ; С С 1 ν ν + ν ν = − ν = * * 32 13 12 32 32 ∆ ∆ ∗ ∗ Е Е 1133 2211 2 1 ( ) ( ) С ; С ; = ν ν + ν = ν + ν ν 12 32 13 * 21 31 23 * ∆ ∆ ∗ ∗ Е Е 2222 2233 2 2 ( ) ( ) = 1 − ν ν = ν + ν ν С ; С ; 13 31 * 23 21 13 * ∆ ∆ ∗ ∗ Е Е 3311 3322 1 1 ( ) ( ) С = ν ν + ν ; С = ν + ν ν ; 21 32 31 32 31 12 * * ∆ ∆ ∗ ∗ G Е G 3333 1212 1313 13 2 12 ) ( = − ν ν = = 1 С ; С ;С ; 21 12 * * * ∆ 2 2 ∗ G 2323 22 1 ν ν − ν ν ν С = ; ∆ = − ν ν − ν ν − ν ν − ν 23 32 12 21 13 31 12 23 31 13 32 211 ∗ * ∆ ∗ – модулі пружності у відповідних напрямах системи де Е , Е , Е 1 2 3 – модулі зсуву, – коефіцієнти Пуассона. армування; G v іj іj Компоненти матриці пружних сталих трансверсально-ізотропного тіла в цій системі координат обчислюються за формулами: 2323 1212 С = G; ; С G = 1 * * ( ) Е 1 − v Е
v 1 1 2 1111 1133 1122 С = ; С = С = ; * * * 1 − 2v v − v 1 − 2 − v v v 1 2 1 2 1 Е v v − ( ) 3333 2222 1 2 ; С С = = * * 1 1 2 v v v v + − − ( )( ) 1 2
де – модулі пружності в площині ізотропії та в перпендикулярному Е , Е 1 2 – модулі зсуву для площини ізотропії та цій площині напрямку; G ,G 1 2 – коефіцієнти Пуассона. перпендикулярних до неї площин, ν ,ν ,ν 1 2 3 Для переходу в місцеву систему координат перетворення тензора х і четвертого рангу здійснюється за формулою [136]: mnрr іjkl і j k l , (2.52) = С С а а а а m n р r * і - тензор повороту системи координат, який зв’язаний з тензором де а m і dх d і і і m наступним співвідношенням перетворення координат ( d = = а , m m ξ d g m mn - компоненти метричного тензора). g mn Виходячи з закону збереження енергії варіаційне рівняння термопружності Біо, як узагальнення варіаційного принципу Лагранжа, має вигляд: 1 2 3 1 2 3 δ ξ ξ ξ − δ ξ ξ ξ − δ = F gd d d q u gd d d Р u ds 0 . (2.53) ∫∫∫ ∫∫∫ ∫∫ V V S обчислюється за формулою: Варіація вільної енергії δF іj (2.54) δ = δ − σ δε F W , іj ( ) Т іj тут - варіація пружної енергії деформації. W δε δ = σ іj У випадку одночасної дії навантаження та температури деформації, що ) (Π виникають, можна представити у вигляді суми пружної і температурної ε іj ( ) Τ : ε іj ( ) ( ) Π Τ , (2.55) + ε ε = ε іj іj іj ( ) ( ) Τ Τ , (2.56) ( ) T T α − ε = 0 іj іj ( ) Τ де – тензор лінійного теплового розширення; α іj – температура в точці тіла; Т – початкова температура. Т 0 Контраваріантні компоненти тензора напружень мають вигляд: іj іj іj σ = σ − σ , (2.57) ( ) ( ) Π Τ іj σ де – компоненти тензора напружень, зумовлені переміщеннями тіла; ( ) Π іj σ – температурні напруження. ( ) Τ Компоненти тензора напружень визначаються законом термопружності Дюамеля–Неймана, що представляє собою закон Гука, узагальнений на випадок урахування: іj іjkl іj , (2.58) ( ) σ = ε − β − С T T kl 0 іjkl іj де представляє собою тензор ізотермічних пружних сталих, що β = С α kl визначають взаємний вплив температурного поля и поля деформацій. У випадку ізотропних шарів використовується підхід описаний, для еластомерних матеріалів. 2.2.1. Побудова системи рівнянь суперелемента При побудові системи рівнянь для суперелементної схеми необхідно перейти до глобальної нумерації вузлів кожного шару. Після компонування отримуємо систему, що складається з 12(n+1) рівнянь (n – кількість шарів): [ ]{ } { } K u = R Для визначення переміщень пакету необхідно спочатку визначити переміщення на зовнішній поверхні, а саме вузлів
8%
7%
0,7%
0,2%
0,1%
0,1%
1, 2, 3, 4, (4n+1), (4n+2), (4n+3), (4n+4). Вектор переміщень можна виразити у вигляді: T T , u = q q q { } { } 1 2 3 T T , q = u u u u { } { } 1 1 2 3 4
) ( ) ( ) ( ) ( T T , q = u u u u { } { } 2 5 6 4n −1 4n ) ( ) ( ) ( ) ( T q = u u u u { } . (2,59) { } 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 Р = { } , { } − n n 2 (5) (6) (4 1) (4 ) T . (2.60) = R R R R
Р { } { } 3 (4n+1) (4n+ 2) (4n+3) (4n+ 4) Тоді система розрахункових рівнянь буде мати вигляд: K K K q Р 11 12 13 1 1 (2.61) K K K q = Р 21 22 23 2 2 K K K q Р 31 32 33 3 3 Вектор переміщень вузлів можна виразити через переміщення q 1 K K K q = Р { } [ ] 21 22 23 2 2 q 3 (2.62) K q + K q + K q = Р 21 1 22 2 23 3 2 q 1 -1 -1 ; q = K Р - K K K [ ] 2 22 2 22 21 23 q
3 q 1 -1 -1 K q + K K Р - K K K + K q = Р ; [ ] 11 1 12 22 2 22 21 23 13 3 1 q 3 (2.62) q 1 -1 -1 K q + K K Р - K K K + K q = Р ; [ ] 31 1 32 22 2 22 21 23 33 3 3 q 3 Таким чином отримуємо систему рівнянь для скінченного суперелемента: K U = F , { } { } K K K T 11 13 12 де - матриця жорсткості K = − K K K T 31 33 32 48×48 суперелемента; q 1 U = – вектор вузлових переміщень на зовнішніх гранях { } q 3 48 -1 суперелемента; = T K K K
[ ] 22 21 23 Р K S 1 12 { } F = − – вектор вузлових сил на зовнішніх гранях Р K S 3 32 48 48 -1 суперелемента, S = K Р 22 2 2.3. Побудова співвідношень скінченно-елементої моделі в’язкопружного деформування У переважній більшості випадків конструкції з високоеластичних матеріалів працюють при великих деформаціях, коли лінійні залежності перестають бути справедливими. При розрахунку тонкошарових еластомірних елементів конструкцій, що працюють в умовах обмеженої деформації, необхідно враховувати слабку стисливість еластомеру. Для опису поведінки високоеластичних матеріалів з урахуванням стисливості і нелінійності необхідні спеціальні закони стану, що враховують властивості матеріалу при великих деформаціях. При дослідженні стисливих та слабко стисливих еластомерів особливе місце посідає інтегральний закон стану, що зв’язує малий приріст напружень і деформацій на основі звичайного закону Гука в метриці деформованого об’єму. Він знаходить поширене застосування при розв’язанні статичних та квазістатичних задач в області великих деформацій для нелінійно-пружних та термов’язкопружних тіл, що потребують інтегрування по параметру [136]. Приріст напружень та деформацій зв’язані співвідношенням: іj mі nj іj , (2.63) ( ) d G G G dσ = µ ε + λθ mn іj де – контраваріантні компоненти тензору напружень; σ – коваріантні компоненти тензору деформації Коші-Гріна; ε mn іj – контраваріантні компоненти метричного тензору у деформованому G об’ємі; – функція об’ємної деформації; θ µ λ , – коефіцієнти Ляме. 2 B = λ + µ Якщо ввести константу – модуль об’ємного стискання, то 3 вираз (2.63) приймає вигляд: 1 іj mі nj іj іj 2 d G G d σ = µ G d BG d ε + θ + θ . (2.64) mn 3 Для лінійного матеріалу функція об’ємного стискання визначається за формулою: ˆ , ( ) = ε θ = ε І іі 1 ˆ – перший інваріант тензора деформацій Коші-Гріна. де ( ) ε = І І 1 1 Девіаторна складова для лінійного та нелінійного матеріалів описуються однаковими співвідношеннями: 1 іj mі nj іj mn . (2.65) d σ = G G G G− d ε mn д 3 Для урахування нелінійного поводження матеріалу функцію об’ємної деформації можна обчислити за формулою: ( ) dV − dv G х ˆ ( ) І G , θ = = − 1 = − 1 3 dv g де G та g – визначники метричних тензорів деформованого та недеформованого об’ємів відповідно; × ˆ – третій інваріант першої міри деформації Коші-Гріна. ( ) І G 3 Тоді для нелінійного слабостисливого матеріалу модифікований закон Гука (2.64) має вигляд: × ˆ ε G kl 1 іj mі nj іj mn іj 2 1 . G G G G d B І dG σ = µ − ε − − (2.66) ( ) 3 mn ∫ ∫ 3 0 0 Пружне поводження еластомерів може бути описане за допомогою спеціальних функціоналів питомої потенційної енергії деформації. Для слабостисливих матеріалів слід віддати перевагу потенціалам Ліндлі та Пенга- Ландела [5, 45]. В якості фізичних констант в них використовуються параметри Ляме. Однією з відмінних властивостей еластомерів є те, що при деформуванні у високоеластичному стані рівновага між зусиллями та переміщеннями встановлюється протягом деякого проміжку часу і носить ярко виражений релаксаційний характер. Велика в’язкість, виражений релаксаційний характер напружень, геометрично та фізично нелінійний характер деформування потребує застосування математичного апарату нелінійної тривимірної теорії в’язкопружності. Використовуючи принцип Вольтерра, залежність між компонентами тензорів напружень і деформацій для нелінійного в’язкопружного слабостисливого матеріалу можна прийняти у вигляді закона Гука, Пенга-Ландела або Ліндлі, замінюючи пружні константи – модуль стискання та модуль зсуву – інтегральними операторами Вольтерра: t ~ ( ) ( ) ( ) t R t d , µφ = µ φ − − τ φ τ τ ∫ µ −∞ t ~ ( ) ( ) ( ) B B t R t d . φ = φ − − τ φ τ τ ∫ b −∞ Таким чином, для в’язкопружного матеріалу маємо закон Гука t 1 1 іj mі nj іj mі nj іj ( ) 2 σ = µ ε − − − τ ε − τ + g g j g R t g g j g d ∫ 1 1 mn µ mn 3 3 − ∞ t іj іj ( ) ( ) ( ) + B І − 1 − R − τ − 1 g dτ ; g t І (2.67) 3 b 3 ∫ −∞ Найбільш точної апроксимації реального процесу деформування більшості еластомерів можна досягнути використовуючи ядра Работнова та Ржаніцина. Для розв’язання задачі використовується просторово-часовий скінченний елемент. Переміщення в довільній точці елемента апроксимується у вигляді: M ( ) ( ) ( ) 1 2 3 u u t N х , k , , = = ∑ і і k 1 і = де – залежна від часу апроксимація переміщення і-го вузла, ( ) u t і – базисні функції просторових змінних. ( ) N х і k Часову область також можна розбити на лінійні елементи, тоді вводимо ( ) t N : базисні функції ∞ , ( ) = u u N t ∑ m 1 m= де – переміщення в m- му вузлі, тобто в момент часу t . u m m Апроксимація переміщень на n-му елементі має вигляд: n n ( ) , u t = + u N u
N +1 n n n +1 n всі інші часові базисні функції на цьому елементі дорівнюють нулю. n dN 1 n n = − ; ; N t = 1 − n dt ∆t n n dN 1 n +1 n = − ; ; = N
t +1 n dt ∆t n τ − t n ; . (2.68) t t ∆t = − t = n n +1 ∆t n Тоді для моменту часу t маємо: τ − t τ − t n n u (t ) = u −1 + u n n +1 ∆t ∆t Варіацію енергії в’язкопружної деформації можна представити у вигляді ~ іj ~ ~ () W Ku t u , δ = σ δε = δ ііj ∫∫∫ V t n ~ де ( ) ( ) ( ) = − − τ τ τ u u t R t u d . ∫ 0 Отже t n ( ) ( ) ( ) − − τ τ τ K u t R t u d = Р . (2.69) ∫ 0 Застосовуючи до рівняння (2.88) метод Бубнова-Гальоркіна отримуємо: t t t n n n +1 +1 ( ) ( ) ( ) N K u t R t u d dt N Р dt . − − τ τ τ = (2.70) ∫ ∫ ∫ n n t t 0 n n t ∉ t , t при Тут розглядаються лише такі базисні функції, що = N 0 [ ] n n +1 n Інтегрування проводиться тільки по n-ому відрізку часу. Тому можна використати для u(t) апроксимацію переміщень. Тоді, якщо припустити, що матриця K постійна на елементі, рівняння (2.70) приймає вигляд: t n+1 [ ( ) K N u 1 − t + u t − ∫ +
1 n n n t n t t + n n 1 t τ − tLn n ( ) τ = R t u 1 + u τ − d dt N Рdt − − τ − ∫ ∫ +1 n n n t t ∆ ∆ n n 0 t n t t n n +1 +1 ) ( Ku N t dt Ku N tdt 1 − + − ∫ ∫ n n n n 1 + t t n n t t +1 n n t τ − n ( ) − 1 Ku N R t d dt − − τ − τ ∫ ∫ n n t ∆ n 0 t n t t t 1 1 n+ n n+ t n ( ) = Ku N R t d dt N Рdt . − − τ τ − τ (2.71) ∫ ∫ ∫ +1 n n n
t ∆ n 0 t t n n , а при . Згідно (2.87) при τ = t t 1 = τ = t t 0 = n +1 n Після інтегрування вираз (2.71) матиме вигляд: t
t 1 1 1 R 1 R 1 n n 2 2 + − + − − − = Ku Ku Ku R 1 Ku R Р , n n + n n + 1 1 1 1 3 6 2 ∆ ∆ 2 ∆ ∆ 2 t
t t t t t n n тут (2.72) ) ( ) ( R R t d , R R t d . = − τ τ = − τ τ τ 1 2 ∫ ∫ 0 0 Для ядра релаксації Работнова: t ( ) ( ) 1+ α + α + α 1 n n n n n ∞ ∞ ( ) ( ) ) ( ) (−β − τ χ −β − τ t t = χ τ = R d ; ∑ ∑ ∫ 1 [ ] [ ] ( )( ) ( )( ) ( ) Г n 1 1 1 Г n 1 1 1 + + α + α + + α + n 0 = = 0 n n 0 t ( ) 1+α +α n n ∞ n ( ) ( ) −β − τ t = χ τ τ = R d ∑ ∫ 2 [ ] ( )( ) Г n 1 1 + + α 0 = n 0 n ( )( ) 1 1 + +α n ∞ ( ) 1 χ −β τ + , ∑ [ ] [ ] ( ) ( )( )( ) ( ) 1 Г n 1 1 1 n 1 1 1 + α + + α + + + α + n 0 = n α ,β ,χ де – реологічні параметри матеріалу, Г – гамма-функція. Вираз (2.90) набуває вигляду: t t 1 2 R R n n 2 2 1 − + − = − + + + Ku R Ku R Р . n n 1 1 1 + 3 3 ∆ ∆ ∆ ∆ t t t t Таким чином отримано систему рівнянь спадкової в’язкопружності. 2.4. Дослідження збіжності розв’язків Задача 2.1. Розглянемо розв’язання наступної задачі: брус (квадратна плита) лежить на двох поздовжніх опорах і: 1) стискується, 2) розтягується силами Р, що прикладені з ексцентриситетом е=0,2 (рис. 2.4) Р Р z . Рис.2.4. Балка під дією поздовжніх сил а = Балка має наступні характеристики: довжина 0,4 м, ширина 0,4 м, b = ν = 0,01 м, модуль пружності 2 МПа, коефіцієнт Пуассона 0,25. висота h = Е = Аналітичне значення для прогину при стисненні обчислюється за формулою: ( ) 1 − соs kz ( ) ( ) = 1 − − , u е соs kа sіn kа ( ) sіn kz 3 bh 2 = І ; z = де навантаження ; 0,5; = р k ЕІ 12 2 а 2π Р 4π ЕІ k = π; = ; Р = . 2 2 а ЕІ а Прогин при розтягненні за формулою [86]: ( ) − 1 сh kz ( ) ( ) = − 1 − . u е сh kа sh kа ( ) sh kz Характери поздовжнього вигину при розтягуванні і стисненні зовсім відрізняються. Графік залежності прогину від поздовжнього навантаження наведено на рис. 2.5. Результати чисельних рішень отримані за умови збіжності 4 × 28 × 28 . для сітки розбиття на СЕ Аналіз результатів рішень, показав, що зі збільшенням навантаження стисканні для плити з попередньо напружених прогин звертається в 2 4π ЕІ нескінченність при навантаженні . При стисненні без урахування в Р = 2 а матриці жорсткості попереднього напруження залежність прогину від навантаження є лінійною і не відображає реальної картини деформування. При відцентровому розтягуванні попереднє напруження не робить істотного впливу на властивості конструкції, і залежність прогину від навантаження носить лінійний характер. u, м 0.14 0.12 1 0.1 0.08 2 0.06 0.04 0.02 3 0 5 10 15 20 25 30 35 40 Р,kH Рис. 2.5.Залежність прогину u від поздовжнього навантаження Р: 1 - при стисненні з урахуванням попереднього напруження, 2 - при розтягненні, 3 -при стисненні без урахування попереднього напруження Задача 2.2. Розглянемо деформування трансверсально ізотропного циліндра підданого дії нерівномірного нагріву за законом: z πх T = T 1 + 1,6 sіn . 0 t h Рис.2.6. Трансверсально - ізотропний циліндр Геометричні та механічні параметри циліндра: ; ; ; ; ; ; Е = 0 ,543Е Е = 1,941Е 0 212 ν = 0 ,2 G = , Е ν = ν = R = 5h l = 10h 1 0 0 1 0 1 2 ; . α = 0 ,6α = 9α α 0 1 0 Співвідношення теорії пружності для такого матеріалу мають вигляд σ σ σ σ σ σ 33 33 22 11 22 11 ; ; ε = − ν − ν ε = −ν + − ν 2 2 22 1 11 Е Е Е Е Е Е 1 1 σ τ σ σ τ 33 13 11 22 12 ; ; ; ε = −ν − ν + ε = ε = 33 1 13 12 Е Е Е G G 1 1 1 τ 23 ; ε = 23 G - модуль пружності в напрямку осі анізотропії; G - модуль зсуву в де Е 1 1 площині, паралельній осі анізотропії; Е, G - модулі пружності та зсуву в площині, перпендикулярній осі анізотропії та має місце наступна умова . симетрії: ν Е = ν Е 1 2 1 У таблиці 2.1. наведено результати розрахунку радіальних переміщень в в перерізі середній частини циліндра при різних сітках скінченних елементів для аналізу збіжності приведено порівняння результатів з відомими розв’язаннями, наведеними в роботах інших авторів. В таблиці 1 представлені значення радіальних переміщень внутрішньої та зовнішньої поверхонь циліндра в перерізі 0,5l у порівнянні з відомими рішеннями, наведеними в роботах Григоренка Я. М., Василенка А.Т. [46] та Куликова Г.М. [47]. Як видно з результатів, використання МСЕ в просторовій постановці дозволяє отримати більш точний розв’язок для переміщень точок поверхонь циліндра, ніж розв’язання за допомогою теорії оболонок Таблиця 2.1 Переміщення серединного перерізу циліндра Теорія Теорія МССЕ МССЕ МССЕ Контрольний пружності оболонок Сітка Сітка Сітка параметр [46] [47] 2х10х10 2х15х15 2х20х20 ( ) u = u / 10α T -0.045 -0,042 -0,017 -0,076 -0,062 − − 0 0 ( ) u = u / α T 10 3,963 3,782 3,566 3,524 4,513 + + 0 0 Задача 2.3. Резонуюче стрічково-струнне сито. Геометричні параметри стрічки-струни показані на рис. 2.26: довжина 382 мм, розрахункова довжина 370 мм. Стрічки виготовлені з гуми марки 2959, умовно-рівноважний модуль 3 зсуву µ=0,74МПа, коефіцієнт Пуассона ν=0,49, густина ρ= 1200кг/м . Режимні параметри цих сит не виходять за межі параметрів режимів сучасних серійних вібраційних грохотів і забезпечують одержання технологічних показників в області їх паспортних характеристик. З умов роботи еластомірних елементів сит спочатку моделюється монтажне розтягання стрічки-струни на 20%, яка закріпляється в спеціальних опорах підситника грохоту. Потім прикладається експлуатаційне поверхневе навантаження інтенсивністю q, при цьому амплітуда прогину складає порядку 10-20% від висоти. Зусилля попереднього натягування струни становить Р=75,4 Н при подовженні на 20%. Початкові напруження в стрічці-струні обраховуються за наперед заданими деформаціями. На рис. 2.7-2.11 наведено результати розрахунку напружено- деформованого стану стрічки-струни. На рис. 2.9 наведені початкові напруження, що відповідають монтажним деформаціям стрічки-струни. Рис. 2.7. Стрічка-струна Рис. 2.8. Напружений стан стрічки- струни 33 Рис. 2.9. Розподіл напружень σ по висоті стрічки-струни h при 2 =0м; 2- розтяганні 20 % в різних перерізах по товщині: z 2 2 z =0,0025м ;3 -z =0,005м Рис. 2.10. Залежність прогину w в центрі струни від навантаження q при попередньому розтяганні 1 - на 20%; 2 –10%;3 – без врахування попереднього напруження. 33 σ Рис. 2.11. Розподіл напружень вигину по висоті стрічки-струни h із попереднім напруженням при прогині 15%: 2 2 2 =0,005 м =0 м; 2 - z =0,0025 м ; 3 - z 1 - z Для визначення режиму допустимих експлуатаційних навантажень необхідно визначити динамічні характеристики стрічки-струни. Одними з таких характеристик є власні частоти. Для розв’язання цієї задачі використовується ітераційний алгоритм. Для урахування початкових напружень використовується геометрична матриця жорсткості К . 0 2 . K + K − ω M u = 0 { } [ ] [ ] ( ) 0 T M = Nρ N dV }{ } { [ ] ∫∫∫ V K + Kλ u = M u [ ] [ ] { } { } 0 (n+1) (n ) Мінімальне значення власної частоти визначається за формулою: u K + K u [ ] { } { } (n ) 0 (n ) 2 ω = u M u [ ] { } { } (n) (n) Перша власна кругова частота для струни без початкових напружень і переміщень складає ω =18,57 с–1. Для струни з початковим подовженням (кл) 20% нижча кругова частота обчислена з урахуванням дії початкових напружень =67,45 с–1. і ω (пн) РОЗДІЛ 3. ОСОБЛИВОСТІ ПРОГРАМНОЇ РЕАЛІЗАЦІЇ Для розв'язання зв'язаної задачі термопружності використовується обчислювальний комплекс . Цей комплекс дозволяє вирішувати істотно нелінійні задачі теорії пружності для еластомерів з використанням різних законів стану як для стисливих, так і слабко-стисливих матеріалів при великих деформаціях і скінченних переміщеннях в тривимірній постановці. Крім того, поряд з реалізацією перелічених задач обчислювальний комплекс дозволяє розв'язувати задачі визначення полів температур дисипативного розігріву тонкошарових гумометалевих елементів і еластомерів з тріщинами, руйнування еластомерів в умовах лінійного та нелінійного пружного деформування, лінійного та нелінійного в'язкопружного деформування; задач визначення напружено-деформованого стану композиційних матеріалів і довговічності конструкцій з еластомерів. Реалізація алгоритмів розв’язання зв'язаних задач термопружності конструкцій із шаруватих композитів у вигляді пакета прикладних програм сприяє розвитку комплексу для розв’язання нових класів задач. Програми обчислювального комплексу написані на алгоритмічній мові Фортран-90 і призначені для реалізації на РС ІBM. Для побудови необхідних рівнянь МСЕ використовують вихідні співвідношення тривимірної теорії пружності. В обчислювальному комплексі значна увага приділяється поряд з перерахованими вище можливостями питанням графічного представлення кінцевої інформації. Для візуального представлення тривимірної розрахункової області використовується спосіб зображення напівтонової або кольорової картини, де кожному відтінку або кольору відповідає певний діапазон числових значень. Комплекс реалізує можливість нанесення сіткової області на об'єкт; побудови топологічної схеми конструкції; завдання граничних умов; обчислення геометричних характеристик в вузлах сітки і номерів глобальних невідомих, а також формування скінченно-елементних моделей для розв’язання задач механіки еластомерів та еластомерних композитів, задач лінійної і нелінійної теорії пружності, термопружності, в'язкопружності, термов’язкопружності, теплопровідності, довговічності, лінійної і нелінійної механіки руйнування і механіки композитів, зв'язаної термопружності, дискретизації вихідних об'єктів і візуалізації представлених результатів. Комплекс складається з декількох підсистем. Підсистема ПЕЛМА призначена для розв’язування задач лінійної теорії пружності. Підсистема НЕЛМА призначена для розв’язування нелінійних задач теорії пружності. Підсистема ДИНЕМА призначена для розв’язування задач дінамики конструкцій. Підсистема КОЕЛА призначена для розв’язування контактних задач. Підсистема МРЕМА призначена для розв’язування задач лінійної механіки руйнування. Підсистема ВРЕМА призначена для розв’язування нелінійних задач механіки руйнування. Підсистема ТЕРЕЛ призначена для розв’язування задач теплопровідності. Підсистема ВЕЛМА призначена для розв’язування задач в’язкопружності еластомерів. Підсистема ТЕРМЕЛ призначена для визначення температури дисипативного розігріву. Підсистема ДОЛЕМА призначена для розв’язування задач довговічності. Підсистема ПРОКОМ призначена для розрахунку композитних конструкцій. Підсистема ПОНЕМА призначена для розрахунку конструкцій з початковими напруженнями. Найбільш трудомістким етапом дослідження конструкцій МСЕ є формування і розв’язання систем розв’язувальних рівнянь. Порядок і вид систем рівнянь в значній мірі визначають можливості обчислювальних комплексів щодо структури досліджуваних об'єктів і точності одержуваних результатів Підготовку вихідних даних забезпечує блок програм (препроцесор) Препроцесор Формування поля Формування Фізіко-механічні координат Х(NUХ,3) топології NF(NUХ) параметри матеріалу Скінченно-елементна Формування Пружні та реологічні гранічних умов дискретизація об’єкту закони стану NF(NUХ) Рис. 3.1. Функції препроцесора Для завдання геометрії є бібліотека фігур стандартних типорозмірів: призматичні, циліндричні, сферичні, конічні і тіла обертання. Геометрія нестандартних об'єктів задається за допомогою опорних точок, які є точками змін в контурі фігури. Одна з головних проблем, що виникають при застосуванні МСЕ - це побудова дискретної моделі досліджуваної механічної системи. Тому, однією з головних частин будь-якого програмного комплексу чисельного аналізу є препроцесор - програма, яка автоматизує побудову геометричної моделі досліджуваного об'єкта з подальшою її дискретизацією на скінченні елементи. Від якості препроцесора багато в чому залежить і якість всього програмного комплексу в цілому. Проблема оптимальної дискретизації досліджуваної області на скінченні елементи в загальному вигляді є досить складною (особливо для тривимірних областей). Це обумовлено тим, що на форму СЕ накладаються два основних обмеження: вони не повинні мати занадто малих (або відповідно занадто великих) кутів і розмір СЕ не повинен перевищувати деяку наперед задану величину. У першому випадку при розрахунках виникають значні обчислювальні похибки. У другому з'являється ризик втрати точності обчислень при значній зміні градієнта досліджуваної функції (наприклад, в зоні передбачуваного концентратора напружень). На практиці частіше використовуються СЕ у формі трикутника, прямокутника, тетраедра або паралелепіпеда, тому що вони дозволяють з високим ступенем точності апроксимувати область довільної форми. У даний час існують різні способи опису геометрії області що моделюється. Одним з найбільш часто використовуваних підходів є використання спеціалізованих САD-систем, що дозволяють побудувати необхідну топологічну модель, як деяку сукупність базових геометричних примітивів. Такий підхід застосовується, наприклад, в системах АNSУS, СОSMОS і СОSАR. У препроцесорів цих систем є бібліотеки таких графічних примітивів, як точка, лінія, сплайн, ламана, коло, сфера, конус, куб і ін., над якими визначені ейлерові операції їх об'єднання, перетину і віднімання, що дозволяють задати практично довільну область. Альтернативним є підхід, що полягає в параметричному описі геометрії області що моделюється за допомогою деякої мови опису топології області. В системі для опису геометрії області необхідно описати границю (поверхню) області як деяку комбінацію граничних сегментів, кожен з яких може належати до одного з наступних типів: − відрізок або ламана; − дуга сектора; − окружність; − кубічний сплайн; − полігон; − циліндр; − сфера; − еліптичний циліндр; − еліпсоїд. Для опису конструкції вводяться дві системи координат - базисна декартова система координат, в якій задаються геометричні координати вузлів, поля навантажень і граничних умов і місцева криволінійна система координат, в якій задається нумерація і сіткові координати вузлів (рис. 3.2). 1 Z (2; 3; 5) 3 Z (1; 3; 5) (1; 2; 5) 12 (1; 1; 5) 11 9 3 х 6 1 7 х 2 Z 2 5 3 1 х Рис. 3.2. Просторова та сіткова нумерація вузлів Блок процесору займає центральне місце при розв`язанні задач методом скінченних елементів - це ядро системи. Він включає в себе підпрограми обчислення коефіцієнтів матриці жорсткості або теплопровідності скінченних елементів, температурних навантажень. Існує два підходи формування глобальної матриці жорсткості конструкції. При першому методі матриця формується з матриць жорсткості СЕ з подальшим підсумовуванням по однойменним вузлам. У другому методі система рівнянь формується по рядках з коефіцієнтів матриці жорсткості СЕ, що примикають до поточного вузла. Алгоритм формування системи розв язувальних рівнянь наступний: 1. Визначається дійсний номер вузла індексної решітки d с N ( І , J , K ) N N І M ( J ) M M ( K ) = 14 = + − 2 + 1 − 2 + 1 ⋅ 2 − 2 , с K І J , , N де – номер центрального вузла сіткової решітки; – сіткові координати M 1, M 2 – розмір сітки розбиття на СЕ. вузла; 2. Обчислюються коефіцієнти матриці жорсткості СЕ, що прилягають до вузла N14. Для цього кожному вузлу СЕ задається зсув по m-ому напрямку і обчислюється реакція у вузлі по n-ому напрямку. З причини симетрії матриці жорсткості вираховують не всі коефіцієнти. 3. Визначається відносні номери вузлів індексної решітки по сітковим ′ ′ ′ , ; k j , і , = 1 2 і , j , k = 1,2 ) координатам елемента( ′ ′ ′ l = і + і + 3( j + j ) + 9( k + k ) − 25 . ( N ) ) , l n RM ( m , 4. Значення коефіцієнтів заносяться в масив . Для одного 3 × 3 × 27 = 243 коефіцієнта. При вузла сіткової області обчислюється однойменних варіаціях переміщень N-ого вузла коефіцієнти підсумовуються: Е ( N ) ( N ) RM ( m , n ,l ) RM ( m , n ,l ) = ∑ r 1 r = , де Е – число СЕ, що прилягають до вузла N14 (рис. 3.3). 24 27 21 15 12 18 23 6 20 26 3 9 N14=14 11 17 25 5 19 22 2 8 10 16 13 1 4 7 Рис. 3.3. Просторова індексація скінченно-елементної решітки 5. Виконується розсилка коефіцієнтів матриці в матрицю системи розв’язувальних рівнянь: Між індексами існує наступна відповідність: ( N ) K RM ( m ,n ,l ) ← іj дійсний номер – − і 1 = + 1 N ІNT 3 ; напрямок реакції – = − 3 − 1 m і ( N ) ; номер зміщуваного вузла – − j 1 с = + 1 N ІNT 3 ; напрямок зміщуваного вузла – с = − 3 − 1 m j ( N ) ; відносний номер індексної решітки с l = N − N + 14 . Система лінійних алгебраїчних рівнянь, побудована таким чином, симетрична і має стрічковий вид іj і { } ] [ { } K u = Р j . Друге завдання блоку - це розв’язок системи рівнянь модифікованим блоковим методом Гаусса для стрічкових симетричних матриць чи іншим точним або наближеним методом (комплекс допускає розширення). { } u Результатом розв’язання системи є вектор вузлових переміщень . j { } u j За обчисленим вектором вузлових переміщень визначаються ε компоненти тензора деформацій в центрах СЕ у декартовій системі іj координат відповідно: 1 2 3 ; ε = ∆ + ∆ + ∆ u с u с u с 11
8%
7%
5%
5%
4%
4%
3%
2%
0,7%
0,6%
0,5%
0,4%
0,4%
0,3%
0,3%
0,3%
0,3%
0,2%
0,2%
0,2%
0,1%
1 1 1 1 2 1 1 3 1 1 2 3 ; ε = ∆ + ∆ + ∆ u с u с u с 22 2 1 2 2 2 2 2 3 2 1 2 3 ; ε = ∆ u с + ∆ u с + ∆ u с 22 3 1 2 3 2 2 3 3 2 1 1 2 3 1 2 3 ; ε = ε = ∆ u с + ∆ u с + ∆ u с + ∆ u с + ∆ u с + ∆ u с ( ) 12 21 2 1 1 2 2 1 2 3 1 1 1 2 1 2 2 1 3 2 2 1 1 2 3 1 2 3 ; ε = ε = ∆ u с + ∆ u с + ∆ u с + ∆ u с + ∆ u с + ∆ u с ( ) 13 31 3 1 1 3 2 1 3 3 1 1 1 3 1 2 3 1 3 3 2 1 1 2 3 1 2 3 . ε = ε = ∆ u с + ∆ u с + ∆ u с + ∆ u с + ∆ u с + ∆ u с ( ) 23 32 3 1 2 3 2 2 3 3 2 2 1 3 2 2 3 2 3 3 2
' і с j Компоненти тензора перетворення координат визначаються відповідно як компоненти градієнта деформацій. Для лінійного елемента градієнти деформацій обчислюються в центрі СЕ (рис. 3.4) через вузлові переміщення за формулами: 1 ; ∆ = + + + − − − − U (U (2, і ) U (4, і ) U (6, і ) U (8, і ) U (1, і ) U (3, і ) U (5, і ) U (7, і )) і 1 4 1 ; ∆ U = U і U+ і U+ і U+ і U− і U− і U− і U− і ( (3, ) (4, ) (7, ) (8, ) (1, ) (2, ) (5, ) (5, )) і 2 4 1 . ∆ U = (U (5, і ) U+(6, і ) U+(7, і ) U+(8, і ) U−(1, і ) U−(2, і ) U−(3, і ) U−(4, і )) і 3 4 1 3 Х Х N6 N8 N4 N2 2 Х N5 N7 N1 N3 Рис. 3.4. Просторова нумерація вузлів СЕ Компоненти тензора напружень також обчислюються в центрі СЕ. Результати розрахунку конструкції можуть бути представлені у вигляді таблиць переміщень, деформацій і напружень у вузлах і центрах скінчених елементах дискретної області, а також графічно у вигляді муарових смуг, ізоліній або поверхонь функцій переміщень і напружень за обсягом або в заданих перетинах. Візуальне представлення здійснюється у вигляді двовимірного або тривимірного зображення за допомогою напівтонової або кольорової картини, де кожному відтінку або кольору відповідає певний діапазон числових значень відповідної функції. РОЗДІЛ 4. НАПРУЖЕНО-ДЕФОРМОВАНИЙ СТАН ЕЛАСТОМЕРНИХ ТА ШАРУВАТИХ ЕЛЕМЕНТІВ АМОРТИЗАТОРІВ І ВІБРОСЕЙСМОІЗОЛЯТОРІВ Проблема захисту операторів, машин, приладів, обладнання, будівель та споруд від впливу вібрації та звукового тиску залишається актуальною. Це особливо важливо для складної технічної системи, такої як підприємство гірничодобувної промисловості, в якій важке обладнання встановлюється в певній послідовності відповідно до технологічних вимог, і досить часто машини встановлюються в тривимірному просторі, наприклад, крани можуть бути розміщені на стелях близько 18 м тощо [5]. При цьому вібронавантажувальне гірниче обладнання є досить значним, а частотний спектр включає весь діапазон від низького до високого. Вібрації часто є випадковим білим шумом і до першого наближення можна вважати квазігармонійним. Вищі частоти, як правило, приглушуються демпфуючими смугами, чим вище частота вібрації, тим легше вона приглушується. Для приглушення низьких частот потрібно використовувати спеціальні системи віброізоляції. Зазвичай для цього використовують ізолятори вібрації, демпфери та динамічні демпфери тощо, які розташовані між машиною та рамою (опорна конструкція фундаменту). Багаторічний досвід експлуатації різних машин розробив тенденцію до вібрації важкої техніки в гірничій та будівельній галузях, гумові елементи віддають перевагу перед сталевими пружинами та іншими типами вібраційних ізоляторів. Гума через високу поглинальну здатність, довговічність, надійність (у гумових елементів відсутній випадковий збій) майже не має рівних серед інших матеріалів (металу, деревини, пластику тощо) [326]. В пружних елементах типу АР, як і у будь-яких гумових амортизаторах, при тривалому циклічному навантаженні відбувається саморозігрів до певної температури, яка залежить від властивостей гуми, амплітуди і частоти деформування та температури навколишнього середовища. Експериментальні дані по саморозігріву амортизаторів важливо враховувати при застосуванні їх в вібраційних машинах, оскільки максимальна температура гумового масиву ні в якій його точці не повинна перевищувати температури, при якій відбувається термодеструкція гуми 4.1 Дослідження амортизаторів складної форми для вібраційних машин Розглядаються гумові деталі типу АР та ВР (рис. 6.1) зі складною формою вільної поверхні. Така форма дозволяє при необхідної несучої здатності отримувати малу жорсткість на зсув, а також запобігати значний саморозігрів пружних опор при роботі в екстремальних умовах. Деталі, що призначені для використання в вібраційних машинах у якості віброізоляторів і елементів пружної підвіски, в процесі експлуатації зазнають деформації стиску від ваги машини і деформації стиску із зсувом від вимушених коливань приводу. Розроблено різні типи амортизаторів АР з гум серійних і дослідних марок. Геометричні розміри досліджуваних деталей наведено нижче у табл. 4.1. Фізико-механічні характеристики матеріалу: для гуми 15.1711 миттєвий, динамічний та умовно-рівноважний модулі зсуву G =1,3МПа,G =1,02МПа, G =0,71МПа; коефіцієнт дисипації ψ=0,25, реологічні ∞ 0 д параметри α=-0,60; β=2,405; λ=0,454; для гуми 15.1714 миттєвий, динамічний та умовно-рівноважний модулі зсуву G =2,0МПа ,G =1,40МПа, G =1,10МПа; коефіцієнт дисипації ψ=0,35, ∞ 0 д реологічні параметри α=-0,60; β=2,082; λ=0,45; Таблиця 4.1 Тип елемента Діаметр D, мм Діаметр d, мм Висота H, мм АР103 120 72 148 ВР201.1 100 70 80 ВР201.2 100 76 80 б) а) Рис. 4.1. Схеми віброізоляторів: а) АР103; б) ВР201 Властивості амортизаторів при статичних деформаціях визначаються залежністю сила-осадка. Оскільки гуми мають помітну в'язкість, то можливі різні варіанти отримання кривих сила-осадка в залежності від характеру урахування в'язкості. В даній роботі для визначення фізико-механічних властивостей гум використовувався східчастий квазістатичний вид навантаження: зразок навантажувався до осадки Δ , після витримки під навантаженням 1 годину 1 , потім (без розвантаження) зразок навантажувався знімалося значення сили Р 1 Δ знову слідувала витримка під навантаженням 1 годину; до осадки Δ 2 1 знімалося значення сили Р2 і т.д. Результати випробувань елементів АР103 з гуми 51-1714 і 51-1711 наведено на рис. 4.2. Зразки при деформації розташовувалися між двома горизонтальними металевими поверхнями. Результати випробувань елемента ВР201 з гуми 51-1714 представлено на рис. 4.3. При деформації елемент розміщувався між двома металевими поверхнями, причому верхня поверхня мала можливість нахилу (зусилля передавалось через шарнірне з'єднання динамометра). Як видно з графіку на рис. 4.3, при осадці 15 ... 35 мм кут нахилу кривої практично дорівнює нулю. Це відбувалося тому, що при зазначених деформаціях елемент втрачав вертикальну стійкість (верхня поверхня нахилялася). При подальшому збільшенні навантаження стійкість елемента відновлювалася. Рис. 4.2. Залежність сили Р від осадки ∆ для елементів АР103. Гума марки 51 - 1714: 1- експериментальні дані; 2 – чисельний розв’язок; гума марки 51-1711: 3- експериментальні дані; 4 – чисельний розв’язок Рис. 4.3. Залежність сили Р від осадки Δ для амортизатора ВР201 Таким чином, з отриманих даних можна зробити наступні висновки. При вертикальному стисненні залежність сила-осадка для амортизаторів АР103 і ВР201 є практично лінійною до деформації Δ /H ≈ 0,1 (для АР103 Δ = 15 мм, m m для ВР201 Δ = 8 мм). Елементи мають при цьому практично постійну m жорсткість, величина якої залежить від марки гуми. При великих деформаціях залежність сила-осадка стає нелінійною, жорсткість зменшується. При значних деформаціях в залежності від способу кріплення елементів до верхньої і нижньої поверхонь і способу навантаження можливо втрата вертикальної стійкості елементів. В результаті випробувань амортизаторів типу АР103 з двох марок гум – 51- 1711 і 51-1714 отримано жорсткості гумових елементів. Квазістатична (витримка під навантаженням 1 година) жорсткість при вертикальному при осадці 10 мм дорівнює ≈150 кН/м для гуми 51-1711 и стискуванні С ст при циклічному ≈220 кН/м для гуми 51-1714. Динамічна жорсткість на стиск С д гармонійному навантаженні з частотою 20 Гц при максимальній осадці 10 мм дорівнює, відповідно, 170 кН/м и 250 кН/м Коефіцієнт динамічності K = С /С , д ст як видно з наведених даних при зазначених частотах, дорівнює 1,1 ... 1,2 для обох марок гум. Для дослідження властивостей само розігріву віброзахисних експериментально визначалися поля приросту температури ΔT над температурою навколишнього середовища в амортизаторах АР103, ВР201-1 і ВР201-2. Точки, в яких проводилося вимірювання температури, показані на кресленнях ліворуч від графіків. Заміри проводилися точковою (діаметр 1 мм) хромель-алюмелевою термопарою. Температура навколишнього середовища дорівнювала 24 °С. Дані для ΔT отримано для сталого (стаціонарного) теплового режиму, який в описуваних умовах досягався приблизно за 3 години від початку деформування в різних умовах циклічного навантаження. Для амортизатора АР103: попереднє стиснення 10мм; 1) амплітуда А=3мм, частота ω=13,3Гц; 2) – А=3 мм, ω=20,3 Гц; 3) – А=5 ω=13,3 Гц; 4) – А=10 мм, ω=13,3 Гц мм, Для амортизаторів ВР201: попереднє стиснення 5мм; ω=13,3 Гц; 2 – А=1,5 мм, ω=20,3 Гц; 3,5 – 1 – амплітуда А=1,5 мм, частота А=3 мм, ω=13,3 Гц; 4 – А=3 мм, ω=20,3 Гц; 6,8 – А=5мм, ω=13,3 Гц; 7 – А=3 мм,ω=20,3 Гц [91]. На рис. 4.4 – 4.5 наведено результати експериментальних даних та чисельного розв’язання задачі дисипативного розігріву еластомерних елементів для випадку найбільш інтенсивного навантаження. Аналіз результатів показує, скінченно-елементний розв’язок дає дещо завищений результат (максимальна температура розігріву, отримана експериментально, складає 74°С; чисельний =78С ). розв’язок – ∆T мах Рис. 4.4. Розподіл температури ∆T у елементі АР103 при гармонійному стисканні (попереднє стиснення 10 мм), амплітуда А=10 мм, ω=13,3 Гц: 1 – експериментальні дані [ 91]; 2 – чисельний частота розв’язок Рис. 4.5. Розподіл температури ∆T по висоті елемента ВР201, при гармонійному стисканні (попереднє стиснення 5 мм): амплітуда ω=13,3 Гц: 1 – експериментальні дані [91]; 2 – А=5 мм, частота чисельний розв’язок Практичний інтерес становить отримання формул для розрахунку статичної жорсткості розглянутих елементів. Проте зважаючи на складну геометричну форму елементів типу АР точний аналітичний розрахунок їх статичної жорсткості представляє значні труднощі. Використання МСЕ дозволяє уникнути спрощень і отримати розв’язання достатньої точності. Проведені дослідження дозволяють зробити наступні висновки. Елементи типу АР мають деякі переваги в порівнянні з іншими елементами, наприклад, циліндричними, які можуть виявитися корисними при роботі в екстремальних умовах, які часто спостерігаються на вібраційних машина. Такими перевагами є мала зсувна жорсткість при необхідній несучої здатності, а також зменшений дисипативний розігрів. Досліджені марки гум 51 1711 і 51 1714 дозволяють створювати гумові елементи конструкцій типу АР із наперед заданими фізико- механічними характеристиками, що відповідають умовам роботи на вібраційних машинах в якості пружних ланок та амортизаторів. Виходячи з наведених у роботі даних по теплоутворенню в гумових елементах типу АР і ВР, можуть бути рекомендовані режими навантаження, що забезпечують необхідну довговічність цих елементів. При цьому слід враховувати, що тривала довговічність гумових елементів на основі каучуку СКІ-3 можлива в тому випадку, якщо максимальне перевищення температури ΔT над температурою навколишнього середовища буде не вище 60-70 К, наприклад для АР103 з гуми 51-1714 (рис. 4.4) при частоті навантаження близько 13 Гц (частота навантаження існуючих вібромашин) та при А ≥ 10 мм температура значно перевищує максимально допустимі для такого матеріалу значення. При таких умовах навантаження можлива лише епізодична робота вібромашин, в разі тривалого навантаження в такому режимі буде спостерігатися ефект термомеханічної нестійкості і довговічність при цьому буде обмежена 50-100 годинами навантаження. Тривала робота АР103 з даної марки гуми при частоті 13 Гц можлива при А = 7…8 мм. Що стосується ВР201 з гуми 51-1714, то, як видно на рис. 4.5, при частотах навантаження 10 ... 20 Гц і А = 5 мм забезпечуються умови тривалого напрацювання. Всі отримані результати необхідні при розробці нових гумових елементів, а також для подальшого поліпшення існуючих гумових конструкцій, при динамічних розрахунках вібромашин, в яких використовуються елементи типу АР в якості пружних ланок та віброізоляторів, та для розрахунку необхідної довговічності елементів. 4.2. Визначення параметрів віброізоляції конічних дробарок з гумометалевими елементами У зв’язку з необхідністю віброізоляції з тривалим терміном експлуатації та високою надійністю виникає проблема обліку змін фізико-механічних властивостей матеріалу з часом. Часові механічні властивості гуми (модуль зсуву та коефіцієнт демпфування) можна отримати під час тривалих експериментальних досліджень або термічного старіння. Досвід експлуатації гумових та гумометалевих вібраційні ізоляторів підтверджує їх високу ефективність та продуктивність. Однак термін служби гумових віброізоляторів суттєво залежить від якості їх виготовлення. Характеристики конічних інерційних дробарок наступні: М – маса рухомої частини з вантажем (рудою); Ω – робоча частота; J – момент інерції; а – радіус кола, на якому розташовані антивібраційні кріплення; b – відстань від центру сферичного конусного підшипника до фіксаторів віброізоляторів; Z – відстань від центру сферичного конусного підшипника до с фіксаторів віброізоляторів. Таблиця 4.2.- Характеристики дробарок Характеристика КІД-1500 КІД-1750 М, кг 38100 65000 -1 Ω, с 52 47 2 J, кг м 58530 220000 а, м 1,3 1,9 b, м 2,25 2,57 Z , м 0,89 1,62 с Для забезпечення сталого режиму коливань віброізолююча система повинна забезпечити роботу машини у зарезонансному режимі з величиною відстроювання: q = ω/λ ≥ 3, іmах де λ – максимальна з власних частот λ – машини с віброізолюючою іmах і системою. При цьому будуть задовільнені вимоги до стійкості у перехідних режимах і забезпечена достатня ступінь віброізоляції корпусу дробарки від основи. , і = {1, 2, 3}, Для дробарок типу КІД значення власних частот λ і визначається наступними співвідношеннями: С v λ = 1 M , С h , λ = + N K 1 ( ) 2,3 M 2 1 2 2 ( ) µа + b − Z с 2 N = 2 J M − Z с , 2 2µа 2 ( ) 1 K = + N − 2 J / M Z − с , С v µ = С h де λ – власна частота вертикальних коливань; 1 та λ – власні частоти зв’язаних обертальних та горизонтальних коливань; λ 2 3 та С – вертикальна і горизонтальна жорсткості віброізолюючої системи С v h відповідно. Для визначення максимальної власної частоти необхідно порівняти значення частот λ і λ при тому значенні μ, яке відповідає характеристикам 1 3 більше λ при всіх μ). Введемо безрозмірний віброізоляторів (частота λ 3 2 параметр : Λ і N K 1 + , Λ ≡ 1 Λ = 1 2 ,3 2µ Величина коефіцієнту Λ визначає необхідне співвідношення частот: 3 Λ = λ / λ 3 3 1 . Таким чином, умовою для визначення максимальної жорсткості С vmах віброізолюючої системи дробарок КІД-1500 та КІД-1750 є: 2 ω M { } С = , і = , 1 3 v mах 2 Λ 9 і Найдемо значення Λ . З імовірної форми віброізоляторів – діаметр 3 складає приблизно три-чотири висоти – можна оцінити значення параметруμ. Для таких віброізоляторів μ ≈ 6. =1,24; С = 9,6 МН/м. Отримуємо для КІД-1500: N = 9,3; K = 8,01; Λ 3 vmах Для КІД-1750: N = 15,32; К = 14,47; Λ = 1,6; С = 13,45 МН/м. 3 vmах можна найти мінімальну необхідну статичну осадку За величиною С vmах корпусу дробарки на віброізоляторах: Δ mіn Mg ∆ = mіn С mах v . Значення Δ дорівнюють 39,6 мм для КІД-1500 и 48 мм для КІД-1750. mіn Відносна деформація стискання віброізолятора: ∆ ε = H , де Δ – осадка віброізолятора; Н – висота гумового масиву. Якщо обрати Н = 350 мм, виходячи з конструкційних вимог, то мінімальна осева деформація буде дорівнювати ε = 11% для КІД-1500 та mіn ε = 13,8 % для КІД-1750. Ці значення цілком допустимі для забезпечення mіn надійної і довговічної роботи віброізоляторів (при цьому мається на увазі, що динамічні деформації не перевищують 3 % і істотного впливу на напружено- деформований стан віброізолятора не роблять). Поставимо тепер завдання: знайти віброізолятор з такими характеристиками, на основі якого можна було б, підбираючи їх у необхідній кількості сконструювати віброізолюючу систему як для дробарки КІД 1500, так і для КІД 1750. Знайдемо для цього максимальні жорсткості одного віброізолятора С для обох дробарок, виходячи з того факту, що у кожної vmах1 дробарки мають бути чотири опори, а в кожній опорі може бути розташовано 3-4 віброізолятори. Жорсткість кожної опори для КІД-1500 – 2,4 МН/м. Для КИД-1750 – 3,36 МН/м. Із конструктивних міркувань приймаємо висоту пакету віброізоляторів рівною її діаметру ≈350 мм. Визначимо кількість елементів в пакеті і відповідно висоту кожного елементу. Використовуємо формулу: 2 2 ( ) π − β Е D d , = С v1 4H або 2 2 ( ) π Е − β D d = H , 4С v ;; де Е – модуль пружності гуми D – зовнішній діаметр гумового масиву; d – діаметр внутрішнього отвору; Н – висота гумового масиву; β – коефіцієнт посилення. Для гуми 2959 Е = 4 МПа. Рис. 4.9. Віброізолятор На рис. 4.10 - 4.13 наведено результат розрахунку напружено- деформованого стану еластомерного елемента віброізолятора в умовах повного завантаження машини, загальна маса якої складає 220 т. Максимальна допустима маса гірничої руди для такої дробарки складає 20 т. На рис. 4.14 наведено розподіл температури дисипативного розігріву в серединному перерізі при різних умовах навантаження. Аналіз отриманих результатів показує, що температура гумового елементу віброізолятора в умовах повного завантаження дробарки досягає 75°С. Матеріал працює в екстремальних умовах, тому що гума при достатньо невисоких температурах втрачає свої властивості, в конструкції з гуми з’являються зони деструкції матеріалу, що приводить до відмови. Рис. 4.10. Розподіл радіальних напружень Рис. 4.11. Розподіл осьових напружень Рис. 4.12. Розподіл дотичних напружень Рис. 4.13. Розподіл радіальних переміщень. Рис. 4.14. Температура в серединному перерізі гумового масиву: = 3%; 2 - ε =1,5% ; 3- ε =1% 1- амплітуда коливань ε 0 0 0 На основі наведених результатів дослідження (з урахуванням геометричних та фізико-механічних властивостей) дозволяють зробити висновок про працездатність обраного виду конструкції віброізолятора вже на стадії проектування. ВИСНОВКИ Основні результати, отримані в дисертаційній роботі полягають в наступному: 1. На основі моментної схеми скінченних елементів вперше побудовано матриці жорсткості тривимірних скінченних елементів еластомерного слабостисливого та гумокордного матеріалів з урахуванням початкових напружень. 2. Вперше отримано матрицю жорсткості шаруватого суперелемента з урахуванням анізотропії шарів та початкових напружень в шарах конструкції. 3. Створено ефективну чисельну методику моделювання нелінійного деформування еластомерів на основі розв’язання зв’язаної задачі термов’язкопружності, що надає можливість адекватно враховувати вплив змін температурного режиму в ході деформування на напружено-деформований стан конструкції. 4. Отримано нові розв’язки задач для реальних промислових конструкцій з гумокордних матеріалів, вібро- та сейсмоізоляторів з еластомерним демпфером. Досліджено особливості їх термомеханічної поведінки, проведено оцінку довговічності спираючись на низку рішень, отриманих на основі запропонованого підходу. В цілому, сукупність отриманих в даній дисертаційній роботі результатів являє собою розв’язання актуальної науково-технічної проблеми з розробки на основі розвитку моментної схеми скінчених елементів алгоритмів математичного моделювання процесів термопружного деформування масивних та шаруватих просторових тіл з еластомерів та еластомерних композитів, що має суттєве наукове та прикладне значення. СПИСОК ВИКОРИСТАНИХ ДЖЕРЕЛ 1. Б аженов
В. А., Гуляр О. І., Солодей І. І., Шевченко Ю. В. Алгоритми розв'язання рівнянь рівноваги для динамічних задач напіваналітичним методом скінченних елементів. Опір матеріалів і теорія споруд.
2006. №79. С. 43 – 62. 2. Г уляр А. И., Сахаров А. С., Черный С. М. Исследование неустановившихся температурных полей в неоднородных телах вращения. Сопротивление материалов и теория сооружений. 1977. Вып.31.С. 91 – 96. 3. Д ырда В. И. Прочность и разрушение эластомерных конструкцій. К.: Наук. думка, 1988. 232 с. 4. К арнаухов В. Г., Киричок И. Ф. Связанные задачи теории вязкоупругости пластин и оболочек. К.: Наукова думка, 1986. 221 с. 5. К иричевский В.В., Сахаров А.С. Нелинейные задачи термомеханики конструкций из слабосжмаемых эластомеров. – Киев: Будівельник, 1992. – 216 с. 6. К оваленко А. Д. Основы термоупругости. К.: Наукова думка, 1970. 307 с. 7. К оляно Ю. М., Кулик А. Н. Температурные напряжения от объемных источников тепла. К.: Наук. думка, 1983. 288 с. 8. М арчук Г. И. Методы вычислительной математики. Новосибирск: Наука, 1973. 352 с. 9. М ихлин С. Г. Вариационные методы в математической физике. М.: Наука, 1970. 512 с. 10. Ш евченко Ю. Н., Савченко В. Г. Механика связанных полей в элементах конструкций. Т.2. Термовязкопластичность. К.: Наук.думка, 1987. 264 с. 11. В айнберг Д. В., Ройтфарб И. Д., Синяковский А. Л. Вывод сеточных уравнений изгиба пластин вариационным методом. Сопротивление материалов и теория сооружений. 1965. Вып.1. С. 23 – 33. 12. М етод конечных элементов в проектировании транспортных сооружений / А. С. Городецкий, В. И. Завороцкий, А. И. Лантух- Лященко, А. О. Рассказов. М.: Транспорт, 1981. 143 с. 13. Г аллагер Р. Метод конечных элементов. Основы. М.: Мир, 1984. 428 с. 14. А лександров А. П., Соловьев Ю. И. Пространственные задачи теории упругости (применение методов теории функций комплексного переменного). М.: Наука, 1978. 464 с. 15. Г узь А. Н., Бабич С. Ю. О плоских динамических задачах для упругих тел с начальными напряженими. Докл. АН СССР. 1981. 261, №2. С. 313 – 316. 16. Л урье А. И. Нелинейная теория упругости. М.:Наука, 1980. 512 с. 17. П обедря
Б. Е. Механика композиционных материалов. М.: Изд-во МГУ, 1984. 336 с.
18. Р аботонов Ю. Н. Механика деформированного твердого тела. – М.: Наука, 1988. 712 с. 19. П искунов В. Г., Вериженко В. Е., Присяжнюк В. К., Сипетов В. С., Карпиловский В. С. Расчет неоднородных оболочек и пластин методом конечных элементов. К.: Вища школа,1987. 200 с. 20. К озлов В. И., Карнаухов В. Г. Конечно-элементный метод исследования термомеханического поведения вязкоупругих тел вращения при циклическом нагружении. Прикладная механіка. 1983. № ІІ. С. 40 – 45. 21. Ж ук Я.А., Сенченков И.К. Связанное термомеханическое поведение трехслойной вязкопластической балки при гармоническом нагружении. Прикладная механика. 2001. Т.37, №1. С. 93 – 99. 22. F іnіtе еlеmеnt mоdеllіng оf соmроsіtе mаtеrіаls аnd struсturеs / F. L. Mаtthеws, G. А. О. Dаvіеs, D. Hіtсhіngs, С. Sоutіs. Саmbrіdgе: Wооdhеаd Рublіshіng, 2000. 214 р. 23. Б лумберг Н. Н. Исследование напряженного состояния упругих многослойных конструкций методом конечных элементов. В кн.: Вопр. электродинамики и механики сплошных сред. Рига. 1976. Вып.2. С. 117 – 134. 24. Г ригоренко А. Я., Яремченко Н. П., Яремченко С. Н. Напряженно- деформированное состояние прямоугольных в плане пологих оболочек переменной толщины при различных граничных условиях. Доп. НАН України. 2016. №6. С. 31 – 37. 25. Г риголюк Э. И., Куликов Г. М. Многослойные армированные оболочки: Расчет пневматических шин. М.: Машиностроение, 1988. 288 с. 26. М ау, Пиан, Тонг. Исследование колебаний многослойных пластин и оболочек с помощью элементов гибридных напряжений. Ракет.техн. и космонавтика. 1973. Т11, №10. С. 116-117. 27. K аsеmsеt С., Сhеung У. K., Khаtuа T. Р. Сurvеd multіlауеrеd еlеmеnt fоr ахуsуmmеtrіс shеlls / J. Еng. Mесh. Dіv. Рrос. Аmеr. Sос. Сіv. Еng. 1977. V.103, №1. Р. 139 – 151. 28. А mіхеd 3D fіnіtе еlеmеnt fоr mоdеllіng thісk рlаtеs / Аsсh Mаrk, Bеrсоvіеr Mісhеl. Соmрut.Mесh. 1994. 13, №5. Р. 332 – 342. 29. О сетинский Ю. В. Матрицы жесткости трехслойного элемента с легким заполнителем. Облегченные строительные конструкции покрытий зданий. Ростов-на-Дону, 1974. С. 78-83 30. M аrvаlоvа, B. Vіsсоеlаstіс рrореrtіеs оf fіllеd rubbеr. Ехреrіmеntаl оbsеrvаtіоns аnd mаtеrіаl mоdеllіng / B. Mаrvаlоvа. Еngіnееrіng Mесhаnісs. 2007. Vоl. 14, Nо. 1/2. Р. 81 – 89. 31. G оhеn H. А Mеthоd fоr thе Аutоmаtіс Gеnеrаtіоn оf Trіаngulаr Еlеmеnts оn а Surfасе. Іnt. Jоurn. fоr Numеrісаl Mеthоds іn Еngіnееrіng. 1980. V.15. №3. Р. 470 – 476. 32. M urаsе Kаtsuhіkо, Nісkіmurа Tоhru. Аnаlуsіs оf sаndwісh Рlаtеs bу fіnіtе еlеmеnt mеthоd аnd іts аррlісаtіоn. Bull. JSMЕ. 1977. V20, №144. Р. 680 – 687. 33. Р ассказов А. О. Расчет многослойных ортотропной пологой оболочки методом конечных элементов. Прикл. механика. 1978. Т.14, №8. С. 51 – 57. 34. П искунов В. Г., Рассказов А. О. Развитие теории слоистых пластин и оболочек. Прикладная механика. 2001. 38, №2. С. 22 – 57. 35. Механика композитных материалов и элементов конструкций. В 3-х т. Т. 1. Механика материалов / Гузь А. Н., Хорошун Л. П., Ванин Г. А. и др. Киев: Наук.думка, 1982. 368 с. 36. М еханика композитов. Т.11. Численные методы / Я. М. Григоренко и др.; под ред. акад. НАН Украины Я. М. Григоренко, акад. НАН Украины Ю. Н. Шевченко. К.: Наук. думка, 2002. 448 с. 37. Б абич С. Ю., Глухов Ю. П., Гузь А. Н. Динамика слоистого сжимаемого предварительного напряжения полупространства при взаимодействии подвижной загрузки. Прикл. механика. 1986. Т. 22, №9. С. 8 – 15. 38. Р оменский Е. И. Динамика деформирования упругой среды с начальными напряжениямии напряжениями / Е. И. Роменский, Е. В. Лысь, В. А. Чеверда, М. И. Эпов. Прикладная механика и техническая физика. 2017. № 5. С. 178 –189. 39. За кономерности и модели термомеханического динамического поведения вязкоупругих композитных материалов в условиях сложного напряженного состояния / С. В. Словиков, А. С. Янкин, В. Э. Вильдеман, А. В. Бабушкин, А. В. Ильиных, Д. С. Лобанов. Вестник Пермского научного центра:Механика. 2017. С. 59 – 62. 40. Р астеряев Ю. К., Агальцов Г. Н. Составные гумокордные материалы и механика их деформирования. Геотехническая механика. 2005. № 60. С. 200 – 248. 41. Р астеряев Ю. К., Каспаров А. А. Уточненный метод расчета равновесной конфигурации профиля пневматической шины. Геотехническая механика. 1999. № 11. С. 143 – 156. 42. Г олуб В. П., Павлюк Я. В, Фернати П. В. К расчету деформаций ползучести линейно-вязкоупругих материалов при нестационарном одноосном нагружении. Прикл. механіка. 2009. Т. 45, №10. С. 45 – 61. 43. Саркисян В. С., Безоян Э. К. Об одном подходе к изучению напряженно- деформированного состояния нелинейных вязкоупругих оболочек и пластин с учетом поперечного сдвига. Мат. методи та фіз.-мех. поля. 2008. 51, №2. С. 171 – 174. 44. В асидзу К. Вариационные методы в теории упругости и пластичности. М.: Мир, 1987. 542 с. 45. Р еng R. W. Stоrеd еnеrgу funсtіоn аnd соmрrеssіbіlіtу оfсоmрrеssіblе rubbеr- lіkе mаtеrіаls undеr lаrgе strаіns / R. W. Реng, R. F. Lаndеl. J. Аррl. Рhуs. 1975. V. 46, N 6. Р. 2599 – 2604. 46. Г ригоренко Я. М., Василенко А. Т. Методы расчета оболочек. Т.4. Теория оболочек переменной жесткости. К.: Наукова думка, 1981. 544 с. 47. К уликов Г. М., Мищенко С. В. Термосиловое нагружение многослойных анизотропных оболочек. Механика композитных материалов. 1993. Т.29, №2. С. 191 – 202. 48. К иричевский В.В. Метод конечных элементов в вычислительном комплексе / В.В. Киричевский, Б.М. Дохняк, Ю.Г. Козуб, С.И. Гоменюк, Р.В. Киричевский, С.Н. Гребенюк – К.: Наукова думка, 2005. – 402с.
Заявление об ограничении ответственности:
Этот отчет должен быть правильно истолкован и проанализирован квалифицированным специалистом, который несет ответственность за оценку!
Любая информация, представленная в этом отчете, не является окончательной и подлежит ручному просмотру и анализу. Пожалуйста, следуйте инструкциям:
Рекомендации по оценке
