Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 12
Текст из файла (страница 12)
Таким образом, для определения интеграла методом Гаусса нужно вычислить в точках интегрирования и сложить и' значений подынтегральной функции, умноженной на весовые коэффи- циенты. Выражение (4.43) можно непосредственно программиро- вать, но прежде введем некоторые обозначения, Обозначим (4.43) (4.44) (и=1, 2, ...,и!). — Л!„ д ду ВЕ = (4.45) (4.46) 77 1.!ЕТ = с1е1 Я (5!, т1!) Н;Н;.
Тогда выражение для блока матрицы градиентов (4,37) четырехугольного конечного элемента с точностью до множителя совпадает с аналогичным выражением (2.64) для треугольного симплекс- элемента. Из этого следует, что результат матричного перемножения в (4.43) при обозначениях (4.44) также с точностью до множителя совпадает с выражением (2.69). Иначе говоря, произведение В нНВ~а' из (4,43) дает тот же результат с точностью до множителя 1/4 (А<'!)', что и аналогичное произведение в (2.67).
Поскольку для вычисления выражения (2.69) ранее была составлена подпрограмма БТ1РР, то возникает вопрос, нельзя ли воспользоваться ею и в данном случае. Подпрограмма ЯТ1РР годится совершенно без изменений для вычисления в одной точке интегрирования блока матрицы жесткости четырехугольного элемента в соответствии с (4.43) при следующих значениях формальных параметров ВЕ и 1)ЕТ.
Двухмерный массив ВЕ определяется по аналогии с (2.70) и предназначен, как и ранее, для хранения ненулевых элементов матрицы градиентов В!'! Элементы массива ВЕ вычисляются в соответствии с (4.38), (4.39), (4.44). Здесь также используется схема хранения нижней симметричной части матрицы жесткости К(') в виде одномерного массива, которая рассматривалась ранее на примерах симплекс-элементов.
Очевидно, подпрограмма сортировки РГ(М5Е элементов блока Кф в одномерный массив ЬЕ ггл ~тся без изменения и в данном случае. Поскольку подпр тграмма ~ТП Е вычисляет блок матрицы жесткости Кф только в одной точке интегрирования, то сразу после вычисления элементы блока Кф должны быть размещены с накоплением в одномерный массив подпрограммой ЕРМАКЕ. Следовательно, совместная работа двух подпрограмм 5Т)ГГ и РКМЗЕ должна быть окаймлена циклами интегрирования по гауссовым точкам по каждой локальной координате. Далее приведен фрагмент программного модуля для вычисления и формирования нижней симметричной части матрицы жесткости изопараметрпческого четырехугольного конечного элемента первого или второго порядка.
Порядок интерполяционного полинома для конечного элемента определяется параметром ХРЕ, которому соответствует число узловых точек в элементе. Таким образом, при ХРЕ = 4 вычисляется матрица жесткости конечного элемента с четырьмя узлами, а при ЯРЕ = 8— с восьмью. Перед началом вычисления матрицы жесткости каждого конечного элемента должна быть сформирована таблица значений координат и весовых коэффициентов для точек интегрирования Гаусса.
СЛ1.1 ТАВОТЕ(14О) ОО 4ИИ 1О=1,Ь(6 СА1Л. ОА1)ЯЯ(16,Х,Н1) 00 4ИИ )6=-1,КО СА1.1. ОА()85() О,У,Н)) 1 СА1-1. Р)ЧГЫ1(14РЕ,Х,У,Х1,У1,Р)Ч1) 2 СА1Л. РЫВИН!()ЧРЕ,Х,У,Х1,71,Р)Ч1) 3 СЛ1.1. РХЗАС(14РЕ,ХРЕ,УРЕ,О(Ч1,РЯС) СЛ1.1. М1ЯУ (ЫС,МОР,ЛЕТ,1.1,1.2) СА1.1, )Ч)ЧВЕЕ(Ч1)Р,)ЧРЕ,1)Ь(1,К)С) ТЕ= — УАШЕ()ЧРЕ,ТРЕ,ЕЬ)1) С1.=РБАБ е У1)УО(ТЕ)ф1+РБАЯ) Ф (1 — 2 * Р('А8)) СМ=У1ЛЧО(ТЕ)/(2э (1+РЖАВ)) ВЕТ=1)ЕТ э Н1э НЯ АРЕА=АВЕА+ОЕТ Я=(3 ж С1.+2 э СМ) э А1.ГА(ТЕ) э ТЕ РО 42И 11.=1,)ЧРЕ Ы= (11.— 1) э Ы1)Е ЭО 43И 1=1,Ы(н'' 78 430 ГЕ(К+1)=ГЕ(И+1)+1)1" 1(1,1Ь) э Я а 1)ЕТ 1)О 420,11 =1,11. СЛ1.1. ЯТ1ГГ(14ЭГ,1Ь..1 Ь,СЬ,СМ,ОХ1,ЯТ,11ЕТ) САЬЬ ГАМЗЕ(МЛГ,1Ь,ЗЬ,БТ,ЗЕ) 420 СОМТ!(Ч(1Е 400 СОМТ1%3Е Поскольку для интегрирования по каждой координате используется одномерный вариант метода Гаусса, то упомянутая таблица представляет собой два одномерных массива, в одном из которых хранятся координаты точек интегрирования, а в другом — значения весовых коэффициентов.
Подпрограмма ТАВЬЕ (см. приложение) формирует по заданному значению формального параметра Хб, который соответствует числу точек интегрирования, массивы координат и весовых коэффициентов. Функциональные возможности подпрограммы ТАВг Е ограничены пятью точками интегрирования, что достаточно для большинства практических задач. Очевидно, вызов подпрограммы ТАВЬЕ должен предшествовать началу цикла по конечным элементам.
'4икл по 1б соответствует интегрированию в (4.43) по координате ~. Подпрограмма ОАЬ)88 (см. приложение) при этом выбирает по номеру 16 из массива координат значение с;, а из массива весовых коэффициентов — значение Н;, которые присваиваются соответственно переменным Х и НЬ В цикле по )6, который соответствует интегрированию в (4.43) по координате т), подпрограмма ОАБ83 выбирает по номеру 36 из массивов координат и весовых коэффициентов значения 1); и Н7 и присваивает их переменным Ъ' и НЯ. Далее алгоритм строится таким образом, чтобы сформировать массив ВЕ на основе представления (4 45) и вычислить значение РЕТ по (4.4б) для выбранной в циклах по 16 и,!О точки интегрирования, определив тем самым фактические параметры для подпрограммы ЯТ1гР.
Процесс формирования массива ВЕ функционально разделен на пять этапов, каждому из которых соответствует свой программный модуль. Подпрограмма РМЕК1 (см. приложение) формирует для выбранной точки интегрирования при значениях Х = ~1 и У = Ч; одномерный массив ГМ, компонентами которого являются значения функций формы (4.28) конечного элемента первого порядка при ЫРЕ = 4 или значения функций формы (4.29) — (4.31) конечного элемента второго порядка при ИРЕ = 8. Формальные параметры Х1 и У1 представляют собой одномерные массивы локальных координат узловых точек соответственно по $ и т).
Программа РХРК1 (см. приложение) формирует для той же точки интегрирования двухмерный массив РЫ, столбцами которого яьляются производные функций формы по локальным координатам. Структура массива ):)Х1 определена соотношением (4.40). В зависимости от значения параметра ЯРЕ массив г)Ы1 соответствует конечному элементу первого или второго порядка. 79 Р'о' = ~ ~ ~"Нвос1е(! Л !с$й1, (4.47) — ! — ! где матрица-столбецначальных деформаций е„для случая плоской деформации дается соотношением (2.72). Прилленяя одномерный метод численного интегрирования Гаусса к интегралу (4.47), получим Геок — Е Х В„" Нва слег 3 ~в!!!, Осо! ° с= — ! !'=! Поскольку подынтегральная функция в (4.47) формально совпадает с видом подынтегральной функции в (2.74), то, очевидно, имеет смысл при программировании выражения (4.48) воспользоваться полученными ранее результатами.
Предварительно отметим тот факт, что при решении задач теории термоупругости ос1 Подпрограмма РМ,1АС (см. приложение) формирует в соответствии с представлением (4.40) и на основе сформированного массива РМ матрицу Якоби, которой отвечает формальный параметр Гс.!С. Параметрам ХРЕ и УРЕ соответствуют одномерные массивы значений глобальных координат узловых точек ХС'Л и л!'с'! на основе соотношения (4.36).
Подпрограмма М[МЧ является стандартной подпрограммой обращения матрицы и вычисления ее определителя. Обратная матрица размещается на месте обращаемой матрицы. Следовательно, результатом работы подпрограммы М1МЧ будет вычисление обратной матрицы Якоби, которая размещается в массиве К,1С, и ее определителя с)ес.1, значение которого присваивается переменной 1ЭЕТ.
Порядок матрицы Якоби равен числу степеней свободы в узле, которое, как и прежде, задается параметром КВР. Параметрам 11 и 1 2 соответствуют вспомогательные одномерные массивы целых чисел, размерность которых совпадает с порядком матрицы. Подпрограмма НОВЕЕ (см. приложение) формирует на основе (4.38) двухмерный массив ВЕ, структура которого дана в (4.45), н размещает его на месте массива 1!!41.
Итак, последовательная работа перечисленных пяти программных модулей, начиная от выбора точки интегрирования Гаусса, завершается формированием массива ВЕ. На основе массива ВЕ следующими двумя циклами по узловым точкам конечного элемента при помощи подпрограмм 8Т(ГГ и гКМЯЕ вычисляется и размещается в массиве Ы симметричная часть матрицы жесткости элемента аналогично тому, как это делалось в случае симплекс-элементов.
Вычисление блока вектора узловых сил, эквивалентных действию распределенного в пределах четырехугольного конечного элемента температурного поля„ сводится к определению интеграла в выражении Контуром четырехугольного конечного элемента являются стороны четырехугольника, которые для элемента второго порядка представляют собой кривые линии также второго порядка.
Для задания любой из сторон конечного элемента можно воспользоваться одномерными функциями формы в виде (4.19) или (4.2), При этом можно отказаться от использования двухмерных функций форм в (4.51), определив уравнение кривой в параметрической форме х = ~ Л' ($)х; (4.52) 7 у = ~' ,Мт©у„(4.53) где у — номер узловой точки , на стороне элемента, а Л' одномерные функции формы. Необходимо помнить о том, чтобы направление обхода узлов е~ на стороне элемента в (4.52) и (4.53) совпадало с направлеа1 х нием обхода всех узлов конечРис. 4.3. К определсниго вектора сди- ного элемента.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
