Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 6
Текст из файла (страница 6)
)' в' К(е) . (2.66) (е) Из (2.66) следует, что матрица К(') симметрична и имеет размерность (6Х6). Программировать непосредственно выражение (2.66) нерационально по той причине, что это потребует сравнительно много ненужных перемножений на нуль. В связи с этим подробнее рассмотрим выражение для блока матрицы жесткости, который определяется так: К(е)„=- ~ В(е)'НВ,(е) Ь (., ~ =-;,), ~). (2.67) (е) Подставим в (2.67) выражения (2.64) и (2.65). Ь О 1 О с() „,, Й). 4А ' са Ьв Ка(з =— (е) П осле выполнения матричных преобразований будем иметь ()л ~~ 21л) ЬелЬа + )лс с() )лЬеес() + рс ЬЗ 1 1 Кее() =- ХЬ()с„+ 1лЬ„са (Х + 21л) с„са + 1лЬ„Ьз~4 А(') (2.68) Перепишем элементы блока Кф матрицы жесткости в несколько иной форме: ~ ХЬ„Ь() + рЬ Ьа ХЬ„св + рс„Ьа 1 ~ Хс„Ь()+ )лЬ с() Хс„са+ (лс ся 14А") á 14(Ь„Ь() + с,са) О 1 1 ЗЗ е ш,ор н.
н. Выражение (2.69) наиболее удобно для программирования, но прежде определим для каждого конечного элемента двухмерный массив ВЕ Ь),~ 1 ВЕ=.~ 12А(~) ' 1с; с1 с),~ (2.70) который попросту предназначен для хранения ненулевых элементов матрицы градиентов В(') . Далее представим Кф в форме Ка)) — 8а)1 + Ра~> (2.71) где 3© — первое слагаемое в матричной сумме (2.б9), а Рф— второе слагаемое в той же матричной сумме. Теперь программирование любого элемента матрицы $<" дается следующим операто- аВ ром Фортрана: ЯТ(1Д) = (В Е(1,11) э ВЕ(3,3 1) э С!+ В Е(3,! 1) э ВЕ(1,З 1) э СМ) а 1) ЕТ Здесь ЬТ вЂ” двухмерный массив матрицы За~)); С1.
= Х, СМ = р„ РЕТ =-- А; двухмерный массив ВЕ определен выражением (2.70), 11 и 11 соответствуют я и р, Для того чтобы получить элемент матрицы блока КД, следует добавить к диагональным элемен гам матрицы Заа соответствую(е) щий элемент диагональной матрицы Р(~а. В целом алгоритм вычисления блока матрицы жесткости К('а) может быть представлен подпрограммой на Фортране ЯЗВ!(ОПТ!ХЕ ЬТ!ГГ(КИ",11.,) 1.,С1.,СМ,ВЕ,ЯТ,ЭЕТ) 01МЕХ510)Ч ВЕ()ЧОГ,!),ЬТ(М1ЭГ,)ЧОГ) 1 ЛО 1И ! — 1,)Ч1)Г 2 1)0 !И )=1,ХВГ 1И ЬТ(1,Ю) = И. 1)0 2И 1=-1,Хг)Г ВО зи К=-1,чг)Г ЗИ ЬТ(1,1)=ЯТ(!,1)+ВЕ(К,!1 ) а ВЕ(К31) ж СМ а!)ЕТ 00 2И 5=1„)Ч1)Г БТ(1,3) ЬТ(1,Д)+(ВЕ(1,11.) а ВЕ(Л,Л1.) а С1.+ а В Е(1,11.) а ВЕ(1,З 1.) а СМ) е ЭЕТ 2И СОХТ1(Ч1)Е КЕТ(1К)Ч ЕМ!) В подпрограмме введены новые параметры: ХРà — число степеней свободы в узле, ВТ вЂ” двухмерный массив матрицы блока Каа.
Введение в список формальных параметров подпрограммы 5Т1РГ значений ЯРЕ и РЕТ будет ясно из последующего. Прежде чем приступить к программированию алгоритма вычисления всей матрицы жесткости К!') элемента, отметим два существенных момента. Во-первых, матрица К<') симметрична и, следова- 34 тельно, вычислению подлежит только симметричная часть. Пусть для определенности это будет нижняя часть матрицы ((х ~ )э). Во-вторых, удобно хранить симметричную часть матрицы К(') в виде одномерного массива, например с~Е (рис. 2.9), в котором элементы симметричной части матрицы К(') расположены по строкам.
Таким образом, требуется подпрограмма, которая сортировала бы элементы каждого вычисленного блока К„р на соответству- (е) ющее место в массиве ЯЕ. Следует иметь в виду, что не все элементы 7 Рис. 2.9. Организация хранения матрицы жесткости двух- мерного симплекс-элемента диагональных блоков должны быть отосланы в массив ЯЕ, а только те, которые принадлежат симметричной части матрицы К(') и ее главной диагонали. Заштрихованные элементы матрицы на рис. 2.9 должны быть проигнорированы подпрограммой сортировки. Перечисленные выше функции сортировки выполняет следующая подпрограмма.
ЯТ)ВКОБТ1)ЧЕ Г)!МБЕ(!ЧРГ,11.,Л1,ЯТ,БЕ) !)1МЕ!ЧБ1ОМ ЯТ()ЧРР,)ЧЭР),ЯЕ(1) Т)О 10 1=!,ИВР е)С вЂ” М!)Р 1Р (1!..ЕЯ.,)1.) МС=! 1С= М!) Р:к (11.— 1)+ 1 00 1О Д=1,КС Л О= И!ЗР:к (Л 1.— 1)+Л )ЧО=--1О е (16 1)/2+36 ! И $Е(Х0)=БЕ(М6)+БТ(1,Л) РЕТ!!КХ ЕМ!) В операторе с меткой 1И используется сложение с накоплением результата, необходимость которого будет видна позже, при рассмотрении конечных элементов высокого порядка. Подпрограммы ЬТ1РР и Р)БОМБЕ дают основу для составления программы вычис- 2* 35 ления симметричной части матрицы жесткости К(е) элемента.
Ниже приведен фрагмент такой программы РО 1И 1!.=1,(ЧРЕ РО 1И 31.=1,11. СА1.1 ЯТ1ГГ(МРГ,!1.,,)1,С1.,СМ,ВЕ,ЯТ„РЕТ) СА1 Е ГРАЧЬЕ(КРГ,11,З 1.,ЯТ,ЯЕ) !И СОМТ1МРЕ Здесь МРŠ— число узловых точек элемента. Очевидно, что перед входом в вызываемые подпрограммы должны быть определены все фактические параметры.
В МКЭ решение задачи термоупругости, т. е, определение напряжений от действия неравномерно распределенного температурного поля сводится к решению задачи изотермической теории упругости путем введения начальных деформаций. Так, например, в случае плоской деформации матрица-столбец начальных деформаций ео() для конечного элемента с температурой Т~'~ определяется следующим образом 1 в(е) = 1 (1+ ъ)(~T(е), (2,72) О В; (е) 1 1'(' =- В1 Н 1 (1 -( ~) ОТ(е) ((Ь. (2.73) В' О (е) Рассмотрим подробнее выражение для блока вектора Г(,), который соответствует узловой точке сс (я =- (', у, й). Получим геоа — ' ) Ва Нво й) = (е) Г (е) ' (е) с '' "1 . 21) е" О 1 Х е(, + 2)( О ! (1+к) ВТ(е) йэ=- О са ЬаЛ О О О Р 1 2,~ (е) (е) 2 (1 — 2ъ) са (2.74) Окончательный результат в соотношении (2.74) справедлив при неизменяющейся в пределах элемента температуре.
Зб где ) — коэффициент Пуассона; Π— коэффициент линейного расширения материала. Тогда можно воспользоваться соотношением (2.15) для определения вектора узловых сил, статически эквивалентных действию начальных деформаций, и записать, что Если температура Т(") изменяется в пределах элемента по линейному закону, то для ее задания можно воспользоваться техникой интерполяции посредством функций формы точно так же, как это делается при интерполяции искомых перемещений. Можно записать Т'» =- Е,Т; -[- Е.Т; + Е,Т~ или в матричной форме Т( Т(е) = [Е, Е, Е,) Т; Т (2.75) где в качестве функций формы используются Е-координаты, Т„ (а = (, 7', А) — температура в узловой точке а. После подстановки (2.75) в (2.73) потребуется вычислить интеграл вида Т 7„= [[ [ [) Т, А,' .
Ть Окончательно выражение для блока вектора Г,, узловых сил, (е) обусловленных линейно изменяющимся в пределах элемента температурным полем, будет (е) Е „ (е) ГеО'е — 2 (! 2е>) 07 е> Я (2.77) где 1 Т() = — (Т;+ Т;+Те). Очевидно, что соотношение (2.77) справедливо и в том случае, когда температура в пределах элемента постоянна. Ниже приведен фрагмент программы для вычисления вектора узловых сил Г,, (е) ЛО 10 11=1,МРЕ (ч=(11.— !) э ь)лг 1)О 10 1=1,Ы!)Е 10 ГЕ(!Ч+1)==ЕЕ(К+1)+ВЕ(!,11) е ():е 1)ЕТ Здесь ГŠ— одномерный массив, в котором хранится вектор Г(е», двухмерный массив определен ранее в (2.70), ЭЕТ =- А('), а ска- лярный множитель Я дается равенством Я =, „ОТ".
Е 37 Т; 7е(» — [ [Е1 Ез Ез) ([() Т> (2.76) А( ) Те. Применение формулы интегрирования (2.55) для интеграла (2.76) дает следующий результат: В представленном выше фрагменте программы опять используется суммирование с накоплением результата, необходимость которого будет видна позже. Вектор узловых сил Ро, статически эквивалентный действию (е) объемных сил 6, определяется на основе (2.17).
В случае плоской задачи теории упругости для симплекс-элемента имеем При условии постоянства объемных сил в пределах элемента и на основании формулы интегрирования (2.55) получим для блока Ре вектора РО~) (е) А (е) ~е (2.78) Вектор узловых сил Р('), статически эквивалентный действию а, распределенных по поверхности элемента сил р, определяется на основе (2.16). Предположим для определенности, что поверхностная нагрузка действует на Рис. 2.!О. К определению вектора единичной нормали к контуру двух- грань элемента, для которой 1 — - О мерного конечного элемента (рис.
2.10). Если поверхностная нагрузка равномерая, то с помощью формулы интегрирования (2.55) легко получить следующее выражение для Р„: (е). О О О О О 1.а О О е'-и (2.79) т,, О О Е. а В случае, когда поверхностная нагрузка изменяется по линейному закону, также можно воспользоваться интерполяцион- 38 Еа У( Рн(' = Еа Л'у ((о) = "(е) Еа Л(» Л л" -7., О— О О О Е, О а О Е, ным соотношением, аналогичным (2.75). Однако для конкретного примера на рис, 2.10 следует учесть, что Е, — О, Тогда имеем Ь и О Е, О (2.80) Подставим (2.80) в (2,79) и, применяя формулу для интегрирова- ния (2.56), получим о 0 Е, Е, О ЕЗ О Е, О Е Ж ~р)й = се> Г у(Е) уй Ея 0 Е2Ез 0 Е'„0 Е2ЕЗ 0 ЕЗ 0 у(Е) /й Е2 Р~и 0 Е,Ез 0 2 0 1 0 2 0 1 0 2 0 1 0 у(о 6 (2.81) (2.83) В случае равенства узловых составляющих поверхностной нагрузки легко видеть, что (2.81) переходит в (2.79).
Обычно в приложениях в качестве распределенной поверхностной нагрузки выступает нормальное давление. В этом случае целесообразно иметь алгоритм и подпрограмму вычисления составляющих нагрузки Р„и Р„по осям координат. Это значит, что требуется определить направляющие косинусы нормали к поверхности элемента, на которой задано давление Р. Будем для определенности разыскивать направляющие косинусы внешней по отношению к элементу нормали ва стороне, для которой Е, — — 0 (рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
