Метод конечных элементов в решении задач теплопроводности Текст научной статьи по специальности «Математика»

УДК 62-278:621.184.64

Доцент С.А. Подгорный, докторант З.А. Меретуков, профессор Е.П. Кошевой, профессор B.C. Косачев

(Кубанский гос. технол. ун-т) кафедра автоматизации производственных процессов

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

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

This work offers to supplement polynomial elements used by the step function of Heaviside and rounding functions which allows to simplify and formalizes the record of test piecewise continuous function applying it for continual problems solution by the method of Galerkin.

Ключевые слова: теплообмен, метод конечных элементов, Галеркин.

Уравнение теплопроводности [1] широко применяется для описания процесса теплопе-реноса в телах классических форм (пластина, цилиндр, сфера), на практике встречаются объекты сложной формы, и задача их описания решается методом конечных элементов [2].

В общем случае рассматривается некоторое семейство функций, определяемых конечным числом параметров. Среди таких функций нет точного решения задачи, однако подбором параметров можно попытаться при -ближенно удовлетворить уравнениям задачи и тем самым построить ее приближенное решение. Специфическим в методе конечных элементов является построение семейства функций, определяемых конечным числом параметров [3-6]. Выберем такое семейство функций u(x) при Xmm < x < Xmax. Интервал Xmm...Xmax представляет собой одномерную область существования решения решаемой задачи, который разбивается на конечное число частей (элементов), соединяющихся между собой и с концами интервала в узловых точках (узлах) X;. В пределах каждого элемента задается функция в виде линейного полинома. Она определяется своими значениями u(X;) в узлах и на концах элемента. Учитывая, что в континуальной задаче функция является непрерывной, то ее значения в каждом узле для соседних элементов должны совпадать. Для этого введем функции округления: LxJ - функция пола, которая определяется как наибольшее целое, меньшее или равное x, а именно LxJ = n ^ x - 1 < n < x; м - функция потолка,

© Подгорный С.А., Меретуков З.А., Кошевой Е.П., Косачев B.C., 2013

которая определяется как наименьшее целое, большее или равное х, а именно |_х] = п ^ х < п < х +1. Используя эти функции округления, введем семейство кусочно-линейных непрерывных функций следующего вида:

M(x) = {L°(x - a)J-f®(x -

c+(d - c)-

x - a b - a

(1)

где Ф(х) - функция Хэвисайда, единичная ступенчатая функция, чьё значение равно нулю для отрицательных аргументов и единице для положительных аргументов; а, Ь - интервал окна функции и(х) на котором она отлична от нуля (а < Ь); с, й - параметры уравнения прямой на интервале а, Ь. Пример такой функции представлен на графике (Рисунок 1).

u(x) 2

a1 1 1 b

I 1 1

Рисунок 1 - Одномерная линейная непрерывная функция

кусочно-

4

2

3

4

5

x

Как видно из представленного графика (Рисунок 1) функция и(х) на границах интервала (а, Ь) обращается в ноль. Используем Ошибка! Источник ссылки не найден. для описания пробной функции на произвольном интервале, задаваемом полушириной к = (Ь - а)/2 одномерного единичного конечного элемента относительно узла X,. В этом случае параметры кусочно-линейной непрерывной функции Ошибка! Источник ссылки не найден. соответственно равны а = X, - к, Ь = X, + к, с = 0, ё = 1, а функция приобретает вид:

% (х) =1_Ф[х - X - к)и - Гф(х - X,,)]}

х + к - X, к

+ {ф(х-X)1-[ф[х-{X, + к)]]}. , (2)

к

Функции ф(х) изображаются в виде ломаных и определяются конечным числом параметров - своими узловыми значениями. На графике (Рисунок 2) показана функция такого семейства.

Мх)

X -к 1 X +к

1 1 1

Рисунок 2 - Кусочно-линейная функция Ошибка! Источник ссылки не найден. для решения континуальной одномерной задачи методом конечных элементов

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

ность более точно удовлетворить уравнениям задачи и тем самым приблизиться к искомому решению. Для линейных задач, когда неизвестные функции и операции над ними входят во все соотношения задачи только в первой степени, метод конечных элементов получил достаточно полное математическое обоснование. В дальнейшем используем линейную задачу, решение которой метод конечных элементов сводит к решению систем линейных алгебраических уравнений. Рассмотрим применение метода для описания одномерного температурного профиля, задаваемого в начальный момент времени (т=0 - начальное условие). В этом случае пробная функция, определяемая уравнением Ошибка! Источник ссылки не найден., представлена линейной комбинацией функций Ошибка! Источник ссылки не найден. с коэффициентами и,=и,(0):

<(х)=£и, -р,(х),

(3)

Для того чтобы и(х,) = и, во всех узлах X,, функции ф(х) должны удовлетворять условиям Ф^) = 1 и щ,(X) = 0 для всех узлов Xj при Кроме того, чтобы выполнялись граничные условия первого рода, следует положить и0 = ип+1 = 0. Метод конечных элементов оперирует в качестве ф,(х) кусочно-полиномиальными функциями, отличными от нуля в пределах небольшого числа элементов вблизи узла X,. Именно это делает метод максимально эффективным. Поскольку и(х) по своему физическому смыслу должна быть непрерывной функцией, выберем у,(х) в виде кусочно-линейных функций, отличных от нуля на двух элементах (Рисунок 2). Каждая такая функция ф(х), , = 1, 2, ..., п, равна единице в X, и нулю во всех остальных узлах. При этом набор функций и(х) будет состоять из непрерывных функций, линейных в пределах элементов с изломами в узлах и определяемых своими узловыми значениями и,, , = 1, 2, ..., п. На концах интервала Xmin...Xmax они обращаются в нуль. Каждую из таких функций можно изобразить в виде ломаной линии. Для определения параметров и,, используемых в уравнении Ошибка! Источник ссылки не найден., сформируем систему линейных алгебраических уравнений методом Галеркина, интегрируя произведение пробной функции на семейство кусочно-линейных непрерывных функций по области существования решения:

2

0

2

4

5

х

тах п

| ф] (хи, (х)

Сх = Пу (х) (х ^¡¡х ,(4)

где]=1,2, ...,п; у(х) - функция начального температурного профиля.

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

Х тах ") У _ У

ш,. = ]>, (х(х)Сх = -• тах ^ т'п , (5)

Х тт

3 п +1

где 1=1,2,..., п.

И над и под главной диагональю:

Ш = ш

Х тах 1 X — X

= (х)-рм (х)Сх = -(6)

6 п + 1

где к=1,2,...,п-1; 1=2,3,.,.,п.

Упростим задачу, полагая, что у(х) =1. В этом случае столбец свободных членов р, представляет собой величину, определяемую интегралом:

х тах х _ X

р. = (х) ах = та* ^

п + 1

(7)

Таким образом, уравнение Ошибка! Источник ссылки не найден. принимает матричную форму:

(1)

2 3 1 6 0 0 0

1 /6 2 3 16 0 0

0 1 6 23 1 6 0

0 0 16 2 3 1 6

0 0 0 1 6 2 3

Здесь матрица является трехдиагональной, и решение может быть получено методом Гаусса, матричным методом (умножением на обратную матрицу) или прогонки. Для систем малой размерности метод решения не существенен, но с увеличением точности решения необходимо увеличивать число конечных элементов на области существования решения, а это приводит к увеличению числа неизвестных в уравнении (1). Для определения влияния числа конечных элементов на интервале Хтт... Хтах на точность начальной аппроксимации решали систему (1) матричным методом в инженерной среде МаШСАБ при различных числах п, определяя среднее значение и(х) на этом интервале. Относительная ошибка при различном числе конечных элементов представлена ниже.

X

X

и

и

2

и

3

и

п-1

и

п

X

—е—6)

О 5 10 15 20 25 30

п

Рисунок 3 - Точность аппроксимации начального температурного профиля в зависимости от числа конечных элементов матричным методом

Из представленных данных видно, что 18 конечных элементах на области существо-

наибольшая точность аппроксимации началь- вания решения. При этом относительная

ного температурного профиля достигается при ошибка аппроксимации начального темпера-12

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

25

а ---а - - 1

1

0 4' 1 4 +аг

. =--

4 + а.

(9)

Р0 =т-, А =

...,Рп_1 =:. 2, (10)

6 ~РП

2 1 4 +«0 1 4 + а п _ 2

Эти величины используются для расчета весовых коэффициентов и,:

и = 6 -Р--1, и 1 =

п 4 + «„ V п1

= ап-1 ■ ип + Рп— , и 1 = «0 • и 2 + Po, (11)

Для определения влияния числа конечных элементов на интервале Хт1П...Хтах на точность начальной аппроксимации решали систему Ошибка! Источник ссылки не найден., Ошибка! Источник ссылки не найден. и Ошибка! Источник ссылки не найден. в инженерной среде МаШСАБ при различных числах п, определяя среднее значение и(х) на этом интервале.

у = 45,480х-° = 0.97{ 371

А

—*—'8. Л Л _ " 11 и о !!■!

д

О 5 10 15 20 25 30 35 40

П, %

Рисунок 4 - Точность аппроксимации начального температурного профиля в зависимости от числа конечных элементов методом прогонки

Относительная ошибка при различном числе конечных элементов представлена выше (Рисунок 4). Из представленных на рисунке данных видно, что метод конечных элементов позволяет снизить относительную ошибку начальной аппроксимации практически до значения менее одного процента при использовании 40 и более конечных элементов. Ис-

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

^ (х )= X и, (т><Р < (х ) , (2)

I = 1

с граничными условиями ^(0) = ^(1) = 0. Однако поскольку уравнение содержит вторую производную по координате, а уже первая производная конечного элемента терпит разрывы непрерывности в узлах, воспользуемся следующим приемом. Обозначим Я(х) = и(т)-ф"(х)-и'(т)-ф(х) невязку исходного дифференциального уравнения. Точное решение дает Щх) = 0. Смягчим выполнение этого условия, потребовав, чтобы оно выполнялось только для п функций, которые совпадают с пробными - и(т)-$(х). Такой прием носит название метода Галёркина. Выполним для невязки интегрирование по частям при условии щ(х) = (Р](х) и ф/0) = (р/1) = 0, тогда получим систему первого порядка, как по временной составляющей, так и по координатной:

хт1х

| к (т) • ф"(х) -«,' (т) • (х)] • Ф] (хУх=

хт!п

- {[-и (Т)-&(х)• ф'](х)-«,'(Т)■ %(х)• ф]=^ (13)

хт1п

Теперь уже в задачу входит и', что дает систему линейных алгебраических уравнений относительно и! вида:

- и }1к'(х)-^ (х)}& =

X тп

= «;(*)• Д^К (х)К (14)

структуры, образованной двумя видами интегралов на главной диагонали матрицы:

m, =

)'Р'(Х = п + 1, (15)

где i=1,2,...,n.

И над и под главной диагональю:

mk ,i +1 = mi ,i -1 = \<P'i(x 1 (x )dx = - n-^1,(i6)

Эта матрица симметричная, что характерно для метода Галёркина. Кроме того, произведение отлично от нуля только при j = i, j = i±1, когда соответствующие два элемента, которые несут на себе пробные функции, перекрываются. Следовательно, матрица в данном случае оказывается трехдиагональной, так как интегрирование совершается только на двух соседних элементах. Решение полученной системы алгебраических уравнений позволяет найти выражение для производных временных составляющих и получить приближенное решение в форме численного интегрирования, например методом Эйлера. Таким образом, для континуальной задачи метод конечных элементов осуществляет приближенный переход к дискретной задаче на основе соответствующих кусочно-полиномиальных функций, отличных от нуля на нескольких соседних элементах, содержащих узел Xi.

Решение (в случае Bi = к) может быть представлено следующей системой дифференциальных уравнений первого порядка:

X

Интегралы правой части Ошибка! Источник ссылки не найден. представлены выражениями Ошибка! Источник ссылки не найден. и Ошибка! Источник ссылки не найден., а интегралы левой части в полученной системе линейных алгебраических уравнений образуют квадратную матрицу ленточного типа трехдиагональной

(-1).

п+1 -К) 0 0 0

-К) п+1 -К) 0 0

0 -К) п+1 0

0 0 -К) п+1 -К)

0 0 0 -К) п+1

«1

и2

« X — X _ max min

п+1

«п-1

«п

23 16 0 0 0

16 23 16 0 0

0 16 23 % 0

0 0 16 23 16

0 0 0 16 23

«П-1

Простейшим вариантом численного интегрирования системы дифференциальных 14

уравнений (37) является метод Эйлера, кото- четной схемой:

рый может быть представлен следующей рас-

4.h) «itaJ w+1 -(-/2) 0 0 0 23 /б 000

«2^ ) «3^ ) «2Ы = «зЫ -(-/2) 0 n+1 -(-/2) -(-/2) n+1 0 -(-/2) 0 0 X — X max min n+! /б 2з 0 /б /б 0 0 23 /б 0

«Äh ) 0 0 -(-/2) n+! -(-/2) 00 /б 2з /б

«h) «ih! 0 0 0 -(-/2) n+! 00 0 /б 23

«^J

«2Ы «зЫ

■Ат

(4)

где Ах - шаг интегрирования системы (37) по времени.

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

ЛИТЕРАТУРА

1 Лыков, А. В. Теория теплопроводности [Текст] / А. В. Лыков. - М.: Высшая школа, 1969.

2 Kosachev, V. S. Using rounding function in the problems of finite-element analysis [Text] / V. S. Kosachev, E. P. Koshevoy, S. A. Podgorny // Studies in mathematical science. 2012. - V. 4. -№ 2.

3 Akin, J. E. Finite element analysis with error estimators [Text] / J. E. Akin. - ButterworthHeinemann, 2005. - 512 p.

4 Chen, Z. X. Finite element methods and their applications [Text] / Z. X. Chen. - Springer, 2005. - 424 p.

5 Solin, P. Partial differential equations and the finite element method [Text] / P. Solin. -Wiley-Interscience, 2006. - 504 p.

6 Thomee, V. Galerkin finite element methods for parabolic problems [Text] / V. Thomee. - Springer, 2006. - 459 p.

REFERENCES

1 Lykov, A. V. Theory of Heat Conduction [Text] / A. V. Lykov. - M.: Vysshaya shkola, 1969.

2 Kosachev, V. S. Using rounding function in the problems of finite-element analysis [Text] / V. S. Kosachev, E. P. Koshevoy, S. A. Podgorny // Studies in mathematical science. 2012. - V. 4. -№ 2.

3 Akin, J. E. Finite element analysis with error estimators [Text] / J. E. Akin. - ButterworthHeinemann, 2005. - 512 p.

4 Chen, Z. X. Finite element methods and their applications [Text] / Z. X. Chen. - Springer, 2005. - 424 p.

5 Solin, P. Partial differential equations and the finite element method [Text] / P. Solin. -Wiley-Interscience, 2006. - 504 p.

6 Thomee, V. Galerkin finite element methods for parabolic problems [Text] / V. Thomee. - Springer, 2006. - 459 p.