Метод конечных элементов в задачах теплопроводности (стр. 5 )

В третьей и четвертой главах были рассмотрены первые два этапа применения МКЭ к решению краевых задач – построение сетки из конечных элементов и получение их базисных функций. В данной главе описаны и, по возможности, проиллюстрированы на конкретных примерах математические процедуры, приводящие к получению расчетных, т. е. программируемых соотношений метода.

5.1 Общее решение краевой задачи методом Галеркина

В рамках метода конечных элементов решение уравнения переноса (1.1.14) можно найти с помощью его вариационной версии или версии метода взвешенных невязок, в частности, метода Галеркина.

В первом случае требуется предварительно сконструировать функционал [2,3,, 3, 4], обладающий экстремальными свойствами, что существенно снижает эффективность вариационного подхода. К тому же функционал, адекватный решаемой задаче, может отсутствовать вообще или не минимизировать искомую функцию в узловых точках. Отсутствие жестких правил регламентации оставляет открытым вопрос о правильности построенного функционала вплоть до стадии проверки размерности получаемых на его основе расчетных соотношений метода. Указанные недостатки наглядно проявляются при работе в системе координат с порядком симметрии [1, 14].

Отсутствие процедуры формирования функционала и использование непосредственно дифференциального уравнения делает методы взвешенных невязок более предпочтительными, особенно если учесть возможность получения с их помощью решения, содержащего варьируемый параметр , т. е. найти общее решение уравнения, справедливое в любой системе координат. 0бщее решение задачи получим с помощью метода Галеркина, для чего достаточно, как указывалось в п. 2.2, найти его для отдельного -го элемента (см. (2.2.2)).

Уравнение переноса (1.1.14) и граничное условие (1.2.7) к нему должны быть записаны для произвольного э-го элемента, что легко достигается приписыванием всем входящим в них функциям и параметрам индекса , указывающего номер элемента. Запишем аппроксимирующую функцию элемента в более удобной матричной форме:

, , (5.1.1)

где [] – матричная строка размером , элементами которой являются базисные функции в узлах элемента; { (ζi)}−вектор-столбец размером значений искомой функции в узлах элемента; − число принадлежащих элементу узлов, ξi , ζi – текущие переменные и координаты узлов, соответственно (i=1,2,3). В дальнейших выкладках для краткости записи индекс e опускаем.

Согласно методу Галеркина приближенное решение уравнения переноса (1.1.14) для -го элемента в общем виде будет описываться выражением (2.2.2), в которое внесен дифференциальный оператор (1.1.14):

, (5.1.2)

где – вектор-столбец базисных функций для элемента с узлами.

Соотношение (5.1.2) содержит дифференциальный оператор 2-го порядка, что не позволяет использовать линейные базисные функции, математически описывающие элементы базового каталога. С целью понижения порядка оператора и включения в уравнение граничных условий (1.2.7), выразим вторую производную функции по обобщенной переменной следующим образом:

. (5.1.3)

Первое слагаемое в (5.1.2) на основании (5.1.3) представим в виде:

. (5.1.4)

Воспользовавшись теоремой Остроградского-Гаусса, заменим в (5.1.4) первый интеграл по объему интегралом по поверхности,охватывающей объем элемента:

.

Подиынтегральное выражение в поверхностном интеграле выразим из предварительно умноженного на граничного условия (1.2.7):

.

.

Используя версию МКЭ (5.1.1), можем написать:

; .

Подставляя полученные путем указанных преобразований результаты в уравнение (5.1.2), решение задачи для элемента e запишем в общем виде:

(5.1.5)

, .

Найденное методом Галеркина общее решение (5.1.5) справедливо для элемента любой размерности и геометрии с произвольным r числом узлов. При вариационном же подходе функционал необходимо конструировать для каждой из систем координат отдельно.

Решение для всей области определения получается суммированием по всем E-элементам элементных вкладов: .

5.2 Матричное представление элементного вклада

Из курса матричной алгебры известно, что произведение матрицы размером (вектор-столбец) на матрицу размером (матричная строка) дает матрицу размером . Поэтому подынтегральные выражения в формуле (5.1.5) для вклада элемента с r узлами в решение задачи представляют собой матрицы размером r x r, или – как последние два члена – вектор-столбцы размером . Так как интегралы берутся в определенных размерами элемента пределах, членами проинтегрированных матриц будут – в конечном итоге – числа, независимо от того, в каком виде удается осуществить интегрирование – в аналитическом или в численном.

Каждая матрица имеет свое традиционное наименование в зависимости от физической природы решаемой задачи. Так, в задачах теплопроводности первый член в (5.1.5) называют объемной частью (интегрирование ведется по объему элемента) матрицы теплопроводности , описываемой выражением

(5.2.1):

выражением:

. (5.2.1)

Верхний индекс показывает, по какой переменной осуществляется дифференцирование базисных функций элемента. Ранг матрицы равен , где –число узлов элемента. У элементов базового каталога r = 4÷8.

Вид производных базисных функций, полученный на основании (1.1.13) – – с учетом значений параметров Ляме в различных системах координат (см. (1.1.10)), представлен в таблице 3.

В матричном представлении первый член в (5.1.5) имеет вид:

, . (5.2.2)

В соответствии с физической природой задачи вместо берем температуру .

Число компонент объемной части матрицы теплопроводности (5.2.1) равно

трем (по числу координат) для естественно ограниченных (S=1) ограниченно симмет-

Таблица 3

Производные базисных функций

S	∂/∂li	dV
i=1	i=2	i=3
1 2 3	∂/∂x ∂/∂r ∂/∂r	∂/∂y (1/r) ∂/∂θ (1/rSinβ) ∂/∂θ	∂/∂z ∂/∂z (1/r) ∂/∂β	dxdydz rdrdθdz r2drdθSinβdβ

.

ричных элементов. У полностью симметричных элементов оно сократится до двух – вследствие азимутальной симметрии -йя компонента исчезнет.

Вторую матрицу в (5.1.5), также умножаемую на вектор-столбец , по аналогии называют матрицей теплопроводности элемента, но ее поверхностной частью (интегрирование ведется по поверхности -ой грани):

, . (5.2.3)

Нижний индекс у матричного произведения базисных функций означает, что базисные функции узлов, не принадлежащих -ой поверхности, должны быть заменены нулями согласно их свойству (4.54.1), а объемные (трехмерные) базисные функции принадлежащих грани узлов – модер­низированы преобразованы, т. е. из объемных превращены в двумерные поверхностные базисные функции, так .как. поверхность описывается уравнением (т. е. нормальна орту ). Очевидно, что и для разных поверхностей выражается по-разному.

Проиллюстрируем сказанное на конкретном примере. Пусть -ой поверхностью является нижняя грань второго элемента каталога и ей присвоен номер 1. Она содержит узлы , и, следовательно, базисные функции остальных узлов . Уравнение плоскости, которой принадлежит первая поверхность: . Это означает, что в базисных функциях узлов этой поверхности текущую переменную нужно заменить на их -e координаты, в силу чего на основании (4.3.3) полином Лагранжа станет равным единице. Для базисной функции на этой грани будем иметь:

, , .

Вторая – верхняя грань – идентична первой, но у нее ,

Вторая – верхняя грань – идентична первой, но у нее , , поэтому .

Остальные грани не являются в общем случае координатными плоскостями, и поэтому базисные функции их узлов нельзя преобразовать; кроме того, , что усложняет процедуру интегрирования. Число компонент объемной части матрицы теплопроводности (5.2.1) равно трем.

Базисные функции тетраэдра также, очевидно, не преобазуются в общем случае, но процедура интегрирования упро­щается переходом к естественной системе координат (см. п. 5.3).

В матричном представлении второй член в решении (5.1.5) имеет вид,

аналогичный (5.2.2):

, (5.2.4)

где - число граней элемента, на которых .

Число компонент поверхностной части матрицы теплопроводности в общем случае равно числу n граней элемента, т. е. n = 4÷6. Для удобства математического представления матрицы (5.1.2) и (5.1.4) объединяют (их ранги одинаковы), что дает:

, (5.2.5)

где –- матрица теплопроводности элемента (матрица жесткости в задачах упругости).

Третий член выражения (5.1.5) принято называть матрицей демпфирования в задачах упругости, и матрицей теплоемкости в задачах теплопроводности:

. (5.2.6)

Как видно из (5.21.6), матрица [] всегда симметрична; (η = (cp ρ )= Cv– объемная теплоемкость материала элемента).

Физический смысл этой матрицы в обеих задачах один и тот же, – она демпфирует (уменьшает) изменения определяемой величины (в частности, температуры) и, тем самым, характеризует инерционные свойства материала элемента (следовательно, и объекта в целом).

Остальные члены,, именуемые объемной и поверхностной частями вектора тепловой нагрузки (вектор нагрузки в задачах упругости), запишутся так:

; (5.2.7)

; (5.2.7)

; где . (5.2.8)

Процедура нахождения поверхностных частей матрицы теплопроводности (5.2.3) и вектора тепловой нагрузки (5.2.8) идентичны. Число компонент поверхностной части вектора тепловой нагрузки равно количеству граней элемента, на которых ≠ 0. В общем случае число компонент равно 4÷6.

Для удобства записи оба вектора объединяют в один вектор тепловой нагрузки:

.

Собирая все члены, получим определяющую элемент систему ( по числу

Собирая все члены, получим определяющую элемент систему (по числу узлов элемента) обыкновенных дифференциальных уравнений первого порядка в их матричном представлении:

. (5.2.9)

Переход от частной производной по времени к обыкновенной объясняется тем, что временная производная рассматривается как функция только координат в каждый фиксированный момент времени. Стационарный случай будет описываться системой алгебраических уравнений:

. (5.2.10)

5.3 Формирование глобальных матриц для исследуемой области

Глобальные матрицы получаются суммированием соответ­ствующих матриц элементов по всем элементам области определения задачи:

; ; , (5.3.1)

; ; , (5.3.1)

где – число элементов, геометрически аппроксимирующих изучаемый объект.

Так как область разбита содержит на R узлов, то система уравнений (5.2.9) превращается в систему R уравнений для R неизвестных значений искомой функции – температуры – в этих узлах:

. (5.3.2)

Очевидно, что при , . Таким образом, ранг глобальной матрицы, а, следовательно, и объем памяти для ее записи, определяется не числом элементов, а количеством узлов. При этом всегда R > E. Поэтому к дискретизации области следует подходить рационально: количество узлов должно обеспечить требуемую точность решения и, в то же время, не превысить слишком завысить выделенный пользователю объем памяти ЭВМ. При дискретизации разных по размерам частей конструкции одним из средств удовлетворения этим требованиям является, как показано в п. 3.2, использование элементов разной геометрии.

В методе конечных элементов каждый элемент рассматривается независимо от других. Поэтому Естественно поэтому при записи его матриц естественно оставлять в ней

ней только те члены, которые относятся к данному элементу, а остальные – нулевые – в его матрицы не включать. Благодаря этому размер матриц элемента равен или , и их называет сокращенными. Такая форма записи существенно экономит объем памяти ЭВМ. Процедуру формирования глобальных матриц рангом легко осуществить, расширив ранг элементных сокращенных матриц до ранга , располагая при этом их члены согласно глобальным номерам узлов элемента (т. н. метод прямой жесткости). При таком способе суммирование Е элементных матриц рангом r, т. е. формирование глобальных матриц, становится тривиальной процедурой (типа известной игры “морской бой “).

Однако использование расширенной матрицы элемента становится крайне неэффективным при большом числе элементов: в памяти ЭВМ пришлось бы хранить матриц рангом .

В эффективных вычислительных программах процедура формирования глобальных матриц осуществляется т. н. методом прямой жесткости, использующего сокращенную матрицу элемента. Сущность метода состоит в следующем.

Строкам и столбцам сокращенной матрицы элемента вместо индексов его узлов приписываются зафиксированные в таблице входных данных (см. стр.40) соответствующие им глобальные номера. Номера строк и столбцов сокращенной матрицы элемента указывают теперь "координаты" каждого члена этой матрицы в сетке глобальной матрицы, что и проиллюстрировано ниже.

а) Сокращенная

матрица элемента

б) Переформированная матрица элемента при ее включении в глобальную

Остальные места в глобальной матрице займут, очевидно, члены сокращенных матриц других элементов. Так как элементы имеет общие узлы, то члены матриц с "координатами" этих узлов будут алгебраически суммироваться. Использование ме-

тода прямой жесткости значительно сокращает загрузку памяти ЭВМ. В реальной про-грамме сокращенные матрицы элементов после их вычисления сразумогут сразу заносияться в соответствующие глобальные и из памяти ЭВМ удаляються.

5.4 Стандартизация матриц элементов

Алгоритм вычислительной программы реализации метода конечных элементов

должен содержать проинтегрированные – стандартизованные – матрицы элементов базового каталога.

Интегрирование матриц элементов каталога, представляющих собой координатные ячейки, не вызывает затруднений, так как в этом случае оно ведется по независимым переменным . Благодаря этому интеграл по объему сводится к произведению трех однократных определенных интегралов, каждый из которых берется по одной из независимых переменных.

Для получения конкретного вида интегрируемой матрицы необходимо проделать следующие операции:

•1) найти конкретные выражения базисных функций всех узлов элемента, ис-

пользуя один из способов математического его описания (см. главу 4) в зави-

симости от типа элемента;

2)• пользуясь таблицей производных (cм. стр. 56), определить производные ба-

зисных функций по текущим переменным;

•3) подставить найденные конкретные функции в соответствующие элементы

матрицы и проинтегрировать их;

• 4) свести в матрицу результаты интегрирования ее элементов.

Согласно (5.2.6) матрица теплоемкости любого элемента всегда симметрична и ее элементами, расположенными на главной и над главной диагональю матрицы, будут интегралы:

и , . (5.4.1)

соответственно; – число узлов элемента. Соотношение (5.4.1) может рассматриваться как тестовое для проверки правильности найденной ранее базисной функции.

Компоненты объемной части матрицы теплопроводности (5.2.1) в разных системах координат описываются следующими конкретными выражениями, полученными с помощью таблицы производных:

1) декартова система координат, ; все три компонента определяются интегралом общего вида:

; ; dV=dxdydz; (5.4.2)

2) цилиндрическая система координат, :

; ; (5.4.3)

; ;

3) сферическая система координат, :

; ;

(5.4.4)

; . (5.4.4)

Подставляя в (5.4.2)–(5.4.4) базисные функции и их производные для конкретного элемента каталога, последующим интегрированием находятся стандартизованные – программируемые сокращенные матрицы теплопроводности. Элементами стандартизованных матриц будут либо числа, либо функции координат элемента. В вычислительной программе значения этих координат берутся из таблицы входных данных. В качестве множителей перед стандартизованными матрицами стоят величины теплофизических коэффициентов из физического каталога.

Процедура нахождения стандартизованных матриц поверхностной части матрицы теплопроводности описана в п. 5.1.

Стандартизация объемной части вектора тепловой нагрузки (5.2.7) сводится к интегрированию базисных функций конкретного элемента каталога. Поверхностная его часть (5.2.8) находится так же, как и поверхностная часть матрицы теплопроводности, но при этом интегрируются при этом сами модернизированные поверхностные базисные функции, а не их произведения.

5.5 Естественная система координат

Когда используется В случае произвольнойая глобальнойая системыа координат, значения узловых координат ограничены только границами области интегрирования. Было бы полезным упрощением, если бы экстремальные значения этих координат были равны –-1, 0, и +1. Этого можно достигнуть выбором локальной (местной) системы координат, привязанной к элементу так, чтобы координаты менялись линейно между нормированными узловыми координатами. Система координат такого типа называется системой естественных координат. Преимущество естественных координат в том, что интегрирование по элементу часто может быть проведено в стандартном аналитическом виде [1, 4, 6].

Естественные координаты могут быть, очевидно, одно - двух - и трехмерными.

В одномерном случае переход к естественной системе координат т осуществляется трансляцией начала системы координат в начало отрезка . Тогда в системе имеем:

,

или, отнеся к длине отрезка:

Рис. 5.1

или, отнеся к длине отрезка,

.

Точно так же, выбирая начало в точке , найдем: .

Из выражений для нормированных координат видно, что это знакомые нам полиномы Лагранжа:

; . (5.4.5)

ТТаким образом, в одномерном случае естественные координаты – это лагранжевы, или -координаты, обладающие, как было показано ранее, всеми свойствами базисных функций. Поэтому аАппроксимирующую функцию элемента можно поэтому записать в естественной системе координат так:

,

поскольку и .

Преимущество веденных -координат в том, что интегрирование можно провести аналитически согласно формуле:

, (5.4.6)

где - длина элемента.

где - длина элемента.

Нормированная кКоордината площади в двумерном случае аналогична нормированной координате длины в одномерном.. Для произвольно выбранной точки в трехузельнловом треугольном элементе такая координата определяется делением площади треугольника А1 (см. рис. 5.2) на площадь А всего треугольника: . Аналогично для остальных частей треугольника:

; и .

Рис. 5.2.

Аналогично для остальных частей треугольника:

и .

Совмещая точку св каждыом из узлов, видим, что в узлах -координаты равны 1, и равны 0 на сторонах, противоположных узлу, что и показано на рисунке. 5.2. Кроме того, очевидно выполнение равенства: