Повна версія

Головна arrow Природознавство arrow СЕТОЧНО-ХАРАКТЕРИСТИЧНІ ЧИСЕЛЬНІ МЕТОДИ

  • Увеличить шрифт
  • Уменьшить шрифт


<<   ЗМІСТ   >>

ПРИКЛАДИ ЧИСЕЛЬНОГО ДОСЛІДЖЕННЯ ДИНАМІЧНИХ ЗАДАЧ МДТТ

На основі викладених в розд. 1 алгоритмів було проведено чисельне дослідження ряду завдань, зокрема, стаціонарних і динамічних плоских задач теорії пружності, двовимірних задач з циліндричною симетрією і урахуванням термічних ефектів для пружного середовища, двовимірних динамічних задач про малі і кінцевих ідеально пружно деформаціях [111, 270-273] . Нижче наведені ряд більш конкретних постановок задач і деякі резлультати цих досліджень, що ілюструють можливості явною і явно-неявної різницевих схем, їх зіставлення один з одним, а також з наявними точними рішеннями. Наведені нижче числові дані обезразмерени наступним чином: лінійні розміри віднесені до деякої характерної для завдання величиною х 0 , щільність - до початкової щільності матеріалу р 0, Компоненти тензора напружень - до деякого характерному значенню про 0 , температура Т - до характерної температурі Т 0 , інші величини відповідно:

і т.д. Індексом р відзначені розмірні величини.

1. Як перший приклад розглянемо плоску стаціонарну задачу про вигин опертої по краях пружною ізотропної плити, кінцевої товщини х 0 і довжини 2х, рівномірним навантаженням (рис. 7.1) [111]. Точне рішення цього завдання (див., Наприклад, [274]) має вигляд

Чисельно завдання вирішувалася в області ABCD (0 <? = Х х < х . 0 <т? = = Х 2 <* 0 ) обмеженою площиною симетрії ВС , з відповідними умовами симетрії на ній

і нижньою межею АВ з рівномірно розподіленим стаціонарної навантаженням про 22 = ~ Я> ° i 2 = 0, а також верхньої вільної кордоном CD (про 2 г = = про 12 = 0) і замикає кордоном AD , на якій задавалося точне рішення з ( 2). Використовувалася описана в розд. 1 постановка, при цьому для даної задачі (f = 0, у = 0, Q = 0, / = 0 - пружний ізотропний матеріал з безрозмірними значеннями параметрів Ляме X = 577, д = = 385) маємо:

Параметри обезразмерени відповідно до (1) (р 0 = 7,8 г / см 3 , х 0 = 1 см, До = 2 • 10 8 Н / м *).

Завдання вирішувалася методом встановлення від деякого початкового розподілу шуканих параметрів і (0, Х, х 2 ) до отримання стаціонарного рішення, тобто до моменту часу, коли із заданою точністю рішення на різних часових шарах збігалося у всіх вузлових точках. В якості початкового розподілу задавалося невозмущенное стан (і (0, х 2 ) = 0) і точне рішення (2). Використовувалася явна схема (1.25) і неявна у напрямку х 2 схема (4.4.7), (4.4.8). Варіювалися також сіткові кроки (20 або 40 вузлів у напрямку х г , 5 або 10 вузлів у напрямку х 2 ). Деякі результати розрахунків при х • = 5, q = 0,1 представлені на рис. 1.2-13.

Як видно з даних на рис. 7.1, на якому наведено поведінку про х i в точці В (*! = 0, х 2 = 0) в залежності від часу /, при фіксованій сітці незалежно від способу завдання початкових даних спостерігається вихід на одні і ті ж стаціонарні значення параметрів (криві 1,2 - відповідно '' нульові "і точні початкові значення). При подрібненні сітки спостерігається збіжність до точного рішення (криві 1-4, відповідні сіток 20X5, 20X5, 20X10 і 40X10). Вплив сіткових параметрів на стаціонарні розподілу про х х (0, х 2 ) і Ог 2 (0, х 2 ) у площині симетрії ВС показано на рис. 7.2, 7.3, штриховими лініями - точне рішення. Тут також видно збіжність до точного рішення при подрібненні сітки, зокрема на сітці 40X10 (h x = 0,125, І 2 = 0,1), максимальні відмінності від точного рішення не перевищують 2-3%. Порівняння розрахунків по явною і явно-неявної різницевим схемам показало, що для однакових значень сіткових параметрів х і І 2 ) розрахунки по обидва схемам практично збігаються, при цьому крок інтегрування т в розрахунках по явно-неявної схеми (4.4.7), (4.4.8) вибирався в 1,5-2 рази більше. ніж у відповідному випадку для явної схеми (1.25).

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

Мал. 7.5

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

В якості другого модельного прикладу розглянемо плоску задачу про нестационарном навантаженні частини верхньої межі пружного шару завтовшки Хо / 2, що лежить (без тертя) на абсолютно жорсткій основі [61, 62]. Область інтегрування (? = Х г > 0, * 0 /2 <т? = Х 2 <* о) зображена на рис. 7.4 і обмежена площиною симетрії Про А з граничними умовами (3), верхньою межею АВС, на частині якої АВ (довжиною х 9 ) граничні умови мають вигляд

а інша частина цієї кордону вільна від навантаження, тобто про 22 _ (* »* ь * о) = = про Х2 (t, Х, Хо) = 0. На нижній межі OD задавалися v 2 ( t % X , х 0 /2) - = 0, а, 2 ( t , х і х 0 /2) = 0. Замикаюча область інтегрування площину CD вибиралася таким чином, щоб за розглянутий час обурення від точки В не досягали цієї межі. В якості вихідних умов задавалося недеформоване стан u (0, х 2 , х 2 ) = 0. Визначальні співвідношення, обезразмеріваніе і параметри Ляме задавалися, як в п. 1. Запис граничних умов в формі (1.34), (1.35) для даного випадку має вигляд

Деякі результати розрахунків при р = М = 1, х 0 = 0,175, х 0 = 1 представлені на рис. 7.4-7.6. Наводиться порівняння результатів, отриманих при h 2 = 0,025 і різних значеннях h х = 0,0333, (темні точки), 0,0167 (хрестики і темні трикутники) і 0,00833 (гуртки). Як видно з даних, представлених на рис. 7.5, на якому показано поведінку компоненти тензора напружень про 2 2 і компоненти v 2 вектора швидкості в точці Х = 0, х 2 = 1 в залежності від часу /, на ділянці АВ верхньої межі спостерігаються коливання а 2 2 біля значення про 2 2 = - р = - 1 і швидкості 1> 2 близько нульового значення з періодом, рівним восьмикратний часу пробігу збурень по товщині шару.

На рис. 7.4 для одного з моментів часу t = 0,03 порівнюється розподіл Оц, і х і про 22 на верхній межі ( х 2 = 1), отримане на різних сітках (Л х = 0,0333, 0,0167, 0,00833) . Видно, що поза невеликій околиці точки В, в якій параметри зазнають розрив, обумовлений розривними граничними умовами, є практична збіжність по сітці. В околиці точки В спостерігається концентрація напруг. Для аналогічної стаціонарної задачі (див., Наприклад, [274]) в околиці цієї точки рішення не обмежена: (Оо ? - 1 / х / дс * - х?)

284

Мал. 7.6

Обшая картина напруженого стану в момент t = 0,03 представлена на рис. 7.6. Показані лінії рівних значень: про 22 ~ рис. 7.6, д, o i2 - рис. 7.6,6, про двох - рис. 7.6, в. Характерним, як зазначалося, є концентрація напружень поблизу точки розриву граничних умов, причому зона великих градієнтів напружень досить значна по глибині.

На рис. 7.7 для одного з моментів часу t - 0,0082 показані ізолінії про 2 2 = const, отримані при розрахунку в аналогічній постановці завдання про раптове і далі постійному за часом і вздовж ОА навантаженні частини кордону пружного півпростору нормальним до цієї межі тиском (плоска задача) [64]. Суцільні лінії на рис. 7.7 і криві 2 на наступному рис. 7.8 відповідають результатам розрахунків по недівергентному варіанту гібридної схеми (4.4.12), штрихові лінії на рис. 7.7

і криві 1 на рис. 7.8 - розрахунок по явною схемою першого порядку точності (1.25). На рис. 7.8 для різних моментів часу наведено розподіл про 2 2 уздовж площини симетрії завдання. Штрихові ламані - точне рішення відповідної одновимірної задачі (коли навантажується вся межа півпростору), рішення якої збігається з рішенням двовимірної задачі уздовж ОВ до моменту приходу до площини ОВ збурень від точки А. Видно, що гібридна схема дозволяє більш детально відтворити дане розривне рішення, проте по обчислювальним витратами вона, природно, поступається простою схемою (1.25).

3. Попередні два приклади ставилися до плоским задачам теорії пружності. Нижче розглядається також двовимірна задача, але з осьової симетрією і з урахуванням термічних ефектів (y ^ OQ ^ O), що моделює дію лазерного випромінювання на слабопоглошаюшую пружну ізотропне середовище. Матеріал і раніше передбачається пружним з визначальним співвідношеннями (1.6), (1.7) при / = 0.

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

В роботі [277] в одновимірної постановці розглянуто задачу про поширення сходяться до осі симетрії хвиль розвантаження, що виникають і рухаються в радіальному напрямку, при частковому поглинанні енергії лазерного імпульсу в циліндричної області. У припущенні т < 11 = R / a , R <х 0 вплив лазерного імпульсу на слабопоглощаю- щий пружний шар моделювався завданням деяких початкових деформацій. Тут г - тривалість імпульсу, R - радіус '' лазерного променя ", Xq - товщина пружного шару, а = [(X + 2д) / р 0 ] 1 ' 2 - швидкість поширення акустичних хвиль в середовищі. Рішення відповідних хвильових рівнянь показало, що хвилі розвантаження, рухомі зі зростанням амплітуди від кордону циліндричного променя, сходяться на осі симетрії до моменту t = 1 1 у що призводить до необмеженого росту тиску на осі. Рішення завдання про сходящейся акустичної хвилі в циліндричної геометрії, яке пророкує, зокрема, необмежене зростання напружень на осі, було отримано раніше в роботі [279].

Нижче розглядається чисельне рішення такого завдання з урахуванням кінцевої товщини зразка ( х 0 ~ R) [64]. Використовується циліндрична система координат Ху = г , х 2 = z, х 3 = <р, так що коефіцієнти Ляме в (1.3), (1.5) мають вигляд Я 1 = Я 2 = 1, Я 3 - г . Область інтегрування (? = = > 0, 0 = х 2 Р ис - 7.9), як і в попередніх двох прикладах,

незмінна в часі (малі пружні деформації з урахуванням термічних ефектів, а й - 0, / = 1,2, / = 1, ..., 7, / = 0, у Ф 0, Q Ф 0) і обмежена віссю ЛВ з умовами симетрії (3) на ній, нижньої ЛЯ і верхньої Z? C межами, на яких задавалося відсутність нормальних і дотичних напружень (а 2 2 = 0, <7i2 = 0), і замикає область інтегрування правої кордоном CD, вибирати таким чином, щоб за розглянутий час обурення від поверхні Ху = R не досягали CD.

Відзначимо, що при використанні систем координат, що мають особливості (наприклад, вісь симетрії jci = г = 0 в циліндричних координатах), замість відповідних умов симетрії (3) можна використовувати асимптотичні рівняння, що виходять з (1.2) граничним переходом г -? 0, що також використовувалося при чисельному рішенні деяких з розглянутих * тут завдань. У формі (1.34), (1.35) граничні умови на AD і ВС мають вигляд СЗ= = | 0, 0, 0, 1, 0, 0, 01, зі! = 3? = {0,0,0,0,

1,0, 0 } ЬЧ = 6? = ь'' = 0.

Об'ємний джерело тепловиділення задавався у вигляді Q (t, x it х 2 ) = = Kq (t, Xi, х 2 ), де до = const - коефіцієнт поглинання, q (г, x iy x 2 ) = 0 При Ху > R , t> 0 і при Ху t> т.

При Ху t <т тепловий потік задавався або однорідним по координаті Xi

або з рішення найпростішого одновимірного рівняння переносу випромінювання, тобто

У представлених нижче розрахункових даних безрозмірні константи матеріалу і характерні розмірні параметри вибиралися в вигляді р 0 = = 2,7 г / см 3 , про про = 6,75 - 10 * 1 дн / см 2 , х 0 = 0,5 см, т 0 = 0,284 • 10 5 к, X = = 0,364, д = 0,318,7 = 0,441, до = 0,5, т = 0,01, з = 1.

На рис. 7.9 для різних моментів часу / суцільними кривими показано розподіл компоненти про х х тензора напружень в залежності від координати х х для одновимірної задачі ( R0 , q (t , х х> х 2 ) з (4), q 0 = 5,56). Тут про хх = про х j / a * віднесено до величини напружень в момент закінчення імпульсу / = т:

Для порівняння штриховий лінією нанесено відповідний розподіл про п з роботи [277] в момент t = 0,5. Певна різниця кривих пояснюється відмінним від розглянутого в [277] деформованим станом в момент Г = т. Видно, що розглянутий тут сеточно-характеристичний метод досить точно відтворює характер рішення практично до моменту '' схлопування "/ = / 1 = 1.

У разі х 0 <2 R хвилі розвантаження, які рухаються від верхньої 2 = * про = О і нижньої 2 = 0) меж шару, взаємодіють в певний момент часу / <1 і викликають помітне перерозподіл напружень у порівнянні з одновимірним випадком (R ^ x 0 ). На рис. 7.10 наведено розподіл в різні моменти часу про 22 = про 22 / о * в залежності від

Мал. 7.11

координати х 2 при х х = 0; простежується вплив двовимірних ефектів (R = х 0 , q (t , л: i, x 2 ) з (4), q 0 = 5,56). Зокрема, тут також спостерігається виникнення великих розтягуючих напружень поблизу осі симетрії х х = 0 в момент '' схлопування "f" 1. При неоднорідному тепловиділення вздовж напрямку х 2 , тобто при використанні (5) замість (4), спостерігаються аналогічні розподілу напружень, кілька деформовані в порівнянні з випадком однорідного тепловиділення.

Для одного з моментів часу t = 0,615 (після моменту взаємодії хвиль розвантаження, що рухаються в напрямку х 2 ) на рис. 7.11 показана загальна картина напруженого стану. Наведено ізолінії про 2 2 = const (рис. 7.11, а), про ц = const (рис. 7.11,6) і про 12 = const (рис. 7.11, в). Для цього моменту часу характерним є поява напруг, що розтягують про х j, про 22 , формування і зростання дотичних напружень про х2 , які локалізуються поблизу кордону променя. Видно область залишкових '' початкових "напруг Оц = - 0,00615.

4. Наступний приклад двовимірної динамічної задачі також стосується малих деформацій (?? • = 0 в (1.5)), але, на відміну від розглянутої вище задачі, тут враховувалися пластичні ефекти в рамках моделі Прандтля-Рейсса за умови пластичності Мізеса (1.7 ) (1Ф 0) [270]. Розглядається задача про дію короткочасної розподіленого навантаження на циліндричну оболонку з зовнішнім радіусом R і товщиною х 0 , підкріплену ребрами жорсткості (рис. 7.12). Термічні ефекти не враховувалися = 0, Q = 0). Завдання вирішувалася з використанням полярної системи координат х х = х 2 = г, х 3 = z так, що коефіцієнти Ляме в (1.3), (1.5) мали вигляд Н х = г , Н 2 = Я 3 = 1. Обезразмеріваніе проводилося

згідно (1), при цьому р 0 = 7,8 г / см 3 , про 0 = 2 • 10 9 дн / см 2 , х 0 = 1 см, X = = 578, ц = 385, к = y R = 15.

Область інтегрування (0 < х х <я, R - х 0 ^ х 2 зображена

на рис. 7.12 і обмежена площинами симетрії А АНН' з граничними умовами симетрії (3), верхньою межею оболонки А '# *, нагруженйой тиском р так, що

В якості вихідних умов задавалося недеформоване стан і (0, х х , х 2 ) = 0.

Деякі результати розрахунків для випадку p (t) = р 0 = 5 при 0 < t < t х = = т / 2 = 0,0136; p (t) = 0 при t > наведені на рис. 7.12-7.13. Тут т = [(X + 2м) / р] " | / 2 - час проходження акустичних збурень поперек оболонки.

На рис. 7.12, 7.13 для моменту часу т < t = 0,043 <2т суцільними кривими показані ізолінії про 12 = const (рис. 7.12) і про 22 = const (рис. 7.13). Видно що чергуються зони стискають і розтягують напруги а 22 , обумовлені різним характером відображення хвилі навантаження від вільних і закріплених ділянок нижньої межі оболонки АН (рис. 7.13). Максимальні значення дотичних напружень формуються поблизу точок розриву граничних умов на внутрішньому кордоні оболонки АН (рис. 7.12). Штриховими лініями на рис. 7.13 позначено межі області пластичних деформацій, побудовані з використанням умов (1.7).

 
<<   ЗМІСТ   >>