Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 8
Текст из файла (страница 8)
= а( 1- (Г. — а() 1 ) -~ — (е — е;) 1,~. (2.1 1 1) Таких( образом, радиус-вектор К из (2.109) посредством соотношения (2.111) является функцией двух координат — Е., и Е,, заданных на поверхности А,, т, е. (е) й = — К (Е., Е,;,). Следуя теории поверхностей, для вектора единичной внешней нормали и) к поверхности А) можно написать формулу (е) (2.112) где Р ь, — — дК1(дЕ,„) и представляет вектор касательной к координатной линии 1.„(и =- 3,2). Далее на основании (2.109) и (2.111) имеем К ь„ =- (х; — х;) е( + (у; †- д;) е~ (- (г; — г;) ез', (2.
113) К ь, = (х~ — х,) е) (-(д), — д() е2 )-(г), — г;) ез. ! Из (2.113) видно, что векторы К ь, и К г, представляют собой не что иное как соответствующим образом направленные стороны треугольника А)( ', т. е. (2.114) откуда следует, по определению векторного произведения, что знаменатель в (2.12) представляет удвоенную площадь треугольника Л,' ): ) К, е, Х К, е,! =2А()').
(2.115) ФВ Числитель в формуле (2.112) определится следующим ооразом; е, е,, е, х,г, у,г, г,~, хг, дг, гг, = Ае) + Ве2+ Сез, где А = (д„— у;) (г, — г;) — (д; — д,) (г~ — г;); В = (г„— г;)(х; — х;) — (г, — г;) (х~ — х;); (2. 116) С = (х, — х;) (у, — у;) — (х, — х,) (у, — у;). Направляющие косинусы нормали к поверхности А~)п получим на основе (2.115), (2.116), (2.112) соз(п,, е)) =- п, е) = А72А,", соз (пь ег) --= п~.е =- В~2А~~е); соз(п,, ез) == и) ез = С/2А~~'.
(2.117) Тогда выражение для определения составляющих нормального давления р на поверхности А) элемента будет иметь вид (е) р„А р,, = В ~ ~(е) р, С (2.118) Для того чтобы формула (2.112) определяла внешнюю к элементу нормаль, нужно внимательно заботиться о соблюдении правила обхода узлов тетраэдра, о котором было упомянуто выше, и исключать соответствующую из 1.-координат в преобразованиях типа (2.110). 2.7. СООТНОШЕНИЯ МКЭ ДЛЯ ТОРОИДАЛЬНОГО КОНЕЧНОГО ЭЛЕМЕНТА В ОСЕСИММЕТРИЧНОЙ ТЕОРИИ УПРУГОСТИ Для решения широкого класса осесимметричпых задач вводится тороидальный симплекс-элемент, образованный вращением плоского треугольного симплекс-элемента вокруг оси симметрии (рис.
2,14). В отличие от плоских и трехмерных конечных элементов, взаимодействие которых осуществляется в узловых точках, взаимодействие осесимметричных тороидальных элементов происходит по узловым окружностям. Осесимметричный тороидальный элемент является двухмерным элементом, интерполяционные соотношения для которого строятся в плоскости, проходящей через ось вращения в системе координат е и г. Следовательно, интерполяционные формулы, полученные для двухмерного плоского симплскс-элемента, пригодны и в данном случае, с тем лишь отличием, что компоненты вектора 49 перемещений и и и являются проекциями последнего соответственно на оси г и г (рис. 2.14).
Определим матричный дифференциальный оператор 0 для осесимметричной задачи теории упругости 0 дт д да (2.119) 0 Г д д дт Тогда по аналогии с (2.63) для осесимметРичной задачи теории упругости блок матрицы градиентов В„' (~ = 1, 1, )е) запишется 1е) 0 д дг Ва Е2Уа. 1 т д дг д де Рис. 2.14. Поперечное сечение тороидального симплекс-элемента (2.120) Подстановка в (2.120) выражений для функций форм из (2 33) " дифференцирование приводят к Ь„' 0 1 2А1') В (е) (2.
121) '"а 2А1е) е Матрица Гука в осесимметричной задаче будет иметь вид Л+2р Л Л О Л Л+2р Л О Л+2и О 0 0 0 р (2.122) Опуская здесь рассуждения, аналогичные тем, которые упоминались в и. 2.5 при определении выражения матрицы жесткости плоского элемента, перейдем к рассмотрению блока матрицы жесткости К®~ тороидального симплекс-элемента.
Поскольку соотношения МКЭ для плоского и тороидального симплекс-элементов во оО многом схожи, то попытаемся выделить в матрице жесткости торо идального элемента те из них, которые присущи только осеси™- тричной задаче. Введем обозначения 0 Ь„О 0 с„ 0 0 0 0 1 2А(е) (2. 123) (е) В, (е) ВОа ~а » 0 п„6 0 0 Л 0 0 Л Л Л~ 2Р 0 0 0 0 0 О 0 0 0 р Л-( 2р, Л Л-г 21( 0 0 0 0 НО н, =- (2 124) Тогда Вее -= Вса + В(ае Н = Но + Н1. Матрицы Во'„' и Н, из (2.123) и (2.124) отвечают только плоской задаче, а матрицы Вф и Н, являются тем дополнением, которое присуще только осесимметричной задаче. Имеем для блока ма трицы жесткости тороидального элемента на основе (2.123) и (2,124) следующее выражение: К('() -- ~ (В(е)'+ В,"„') (Н,,— Н,) (В,',',) — В(')) (1 . (2.128) фе) Интегрируя (2.125) по угловой координате и ччитывая только ненулевые матричные произведения, получим Кф = — 2л ~ Вф Н Вф»(1»Ж+2д ~ (В(('„~ Н)В)() + В(еез Н(Воа + ,((е) А(е) +В,'„'Н(Вфр')»(1»(1г=1А-(-1в (а, (" = ( ! ~).
Заметим, что в подынтегральные выражения в (2.126) входит переменная величина», которая определяется »(е) — Е,»(+ Е,»;+ Ез»и. (2.127) Вычисление первого интеграла в предположении зависимост" (2.127) для» дает 1А = — Воа НОВ()в)У(', (2. 128) где $" = 2л»")А", а»(' =- — (»(+») +»),). Результаты матричных преобразований в (2,128) после подстановки соответствующих выражений из (2.123) и (2.124) показываю~ что интеграл 1„с точностью до скалярного множителя соответст 51 вует матрице жесткости блока плоского симплекс-элемента (2.69).
Это обстоятельство должно быть учтено при программировании вычисления блока матриць. жесткости тороидального элемента. Вычисление второго интеграла в (2.126) осложняется тем, что осесимметричная часть матрицы градиентов В(„тороидального (е) элемента зависит от г. Однако на практике при работе с симплекс- элементом величина интеграла 1а вычисляется при постоянном значении радиуса в центре треугольного сечения г — г'). Тогда ~" г ,(е) ),~(е) 0 т '/П т 'ее Л'р т.са (е) е" (2.129) Очевидно, что для вычисления блока матрицы жесткости К„в (е) тороидального симплекс-элемента нельзя воспользоваться готовой подпрограммой ЯТ1ГГ, поскольку в последней не предусмотрено формирование матрицы 1и из (2.129). Необходимые изменения подпрограммы ЯТ1ГГ сводятся к следующему.
Во-первых, дополнительно в качестве формального параметра требуется ввести одномерный массив КХ, элементами которого являются ненулевые элементы третьей строки матрицы градиентов В(') конечного эле- мента и; м„Ф~ 1 КМ = ,(.) ,(.) ,( ) т т т (2. 130) ясВКОБт1)че ят1елр)РГ,11.„)).,с!.,см,Ве,Кь),ят,Рет) Р1МЕМЯОЫ ВЕ(ХРГз !),ЬТ(ХРГ,МРГ),КХ(!) ят0,1)=(к)ч(11.) э (с1-+2 э см) ~ кщ) 1-)+ е С1. е (КК(11.) е ВЕ(1,31)+ КМ(,! 1.):е ВЕ(1,! 1.))) е РЕТ 5Т(1,2)=С1.э ВЕ(2,)1.) е К)Ч(11) е РЕТ ВТ(2, !) — С1.:е ВЕ(2,11.) е К)Ч(!1) ~ РЕТ ЯТ(2,2) = И.
РО !И !=-1,КРЕ РО 2и К= ),(чРГ 2И ЯТ(!,1)=ВТ(1,1)+ВЕ(К,11) э ВЕ(К )1) эСМэ РЕТ РО )И .)=- ),КРЕ ЬТ(1,))=ЯТ(1,))+(ВЕ(1,11) э ВЕ(5,31) е С1+ е ВЕ(5,11-) е ВЕ(1,)1 ) э СМ) э РЕТ 1И СОМТ1)ЧБЕ КЕТ()КХ Е)ЧР Во-вторых, операторы с метками 1, 2, 10, которые соответствуют зачистке массива ЯТ, следует заменить операторами присваивания элементам массива 5Т значений элементов матрицы К~ч~~, Таким образом, подпрограмма вычисления блока матрицы жесткости тороидального элемента с именем ЯТ1ГА может иметь следующий вид.
1 ОТ('). 1 еО (е) (2.131) Тогда для вектора узловых сил, статически эквивалентных дейст- вию начальных деформаций, можно записать (е] ОТ(') сЬ (2.132) в; в' или для блока вектора Ген = Ва 11е0 (10 = ° г г(е) ܄Π— 2А(') =- 2л 2Аьч О с О А(~) а Ср — ОТ(']г й' Ж (1 — 2ч') ()а О Ь + " 2А(') ОТ(е)г (Ь дг. (1 — 2ъ) (2.134) я (е) Подынтегральная функция в (2,134) довольно сложным обра- зом зависит от г, так как температура в элементе обычно является 53 Для хранения элементов матрицы жесткости К('] тороидальпого элемента, как и для плоского симплекс-элемента вводится одномерный массив ЯЕ.
Поскольку размерность и структура матрицы жесткости тороидального элемента в точности соответствуют размерности и структуре матрицы жесткости плоского элемента, то все рассуждения относительно сортировки блока К а матрицы плоского элемента в одномерный массив ЯЕ справедливы и в данном случае (см. рис. 2.9). Таким образом, подпрограмма ГКМЯЕ без изменений пригодна для формирования элементов матрицы жесткости тороидального элемента в виде одномерного массива 5Е. Ясно, что фрагмент программы для вычисления симметричной частицы матрицы жесткости тороидальпого элемента будет соответствовать фрагменту на стр. 36, с тем лишь отличием, что обращение к подпрограмме ЬТ1ГГ должно быть заменено обращением к подпрограмме ЬТ1ГЛ.
Для решения задач осесимметричной теории термоупругости введем матрицу-столбец начальных деформаций конечного элемента с температурой Т(') функцией радиуса. Однако в целях упрощения интегрирования предположим сначала, что подынтегральная функция постоянна и вычислена при значении г = г~, т. е. в центре треугольного се(е) чен ия элемента . Получим оа Л' а 2А(') (2.135) са 2А(') Формула (2.135) легко программируется и ~а основании введенных ранее двухмерного массива ВЕ (2.70) и одномерного массива ЕК (2.!30) можно написать следующий фрагмент программы на Фортране: Т)О 18 И.=),Я'Е Х=(11.— 1) э МРР РЕ(К+1)= РЕ(~х 1)+(ВЕ(1,11)+КХ(!1)) а Яе ВЕТ 10 РЕ(И+2)=РЕ(К-т-2)+ ВЕ(2,11-) -в Я е!)ЕТ 2т(г(')А", а Я дается Здесь МРŠ— число узлов элемента, Е)ЕТ выражением (2.136) В случае точного интегрирования в (2.134) требуется вычислить два следующих интеграла по площади треугольника: .