Вестник ВГУИТ, №2, 2013

УДК 62-278:621.184.64

Доцент С.А. Подгорный, докторант З.А. Меретуков,

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

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

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

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] широко которая определяется как наименьшее целое,

применяется для описания процесса теплопе- большее или равное x, а именно
реноса в телах классических форм (пластина, x = n ⇔ x < n ≤ x +1. Используя эти функции
цилиндр, сфера), на практике встречаются округления, введем семейство кусочно-линейных
объекты сложной формы, и задача их описания непрерывных функций следующего вида:
решается методом конечных элементов [2].
В общем случае рассматривается неко-  x − a
u (x ) = {Φ(x − a ) − Φ(x − b )}⋅ c + (d − c ) ⋅ , (1)
торое семейство функций, определяемых ко-  b − a 
нечным числом параметров. Среди таких
функций нет точного решения задачи, однако где Φ(x) – функция Хэвисайда, единичная сту-
подбором параметров можно попытаться при- пенчатая функция, чьё значение равно нулю
ближенно удовлетворить уравнениям задачи и для отрицательных аргументов и единице для
тем самым построить ее приближенное реше- положительных аргументов; a, b – интервал
ние. Специфическим в методе конечных эле- окна функции u(x) на котором она отлична от
ментов является построение семейства функ- нуля (a ≤ b); c, d – параметры уравнения пря-
ций, определяемых конечным числом пара- мой на интервале a, b. Пример такой функции
метров [3-6]. Выберем такое семейство функ- представлен на графике (рисунок 1).
ций u(x) при X min ≤ x ≤ X max . Интервал
X min …X max представляет собой одномерную a b
область существования решения решаемой
задачи, который разбивается на конечное чис- 4 d

ло частей (элементов), соединяющихся между

собой и с концами интервала в узловых точках
(узлах) X i . В пределах каждого элемента зада-
ется функция в виде линейного полинома. Она u( x) 2

определяется своими значениями u(X i ) в узлах c

и на концах элемента. Учитывая, что в конти-
нуальной задаче функция является непрерыв-
ной, то ее значения в каждом узле для сосед- 0

них элементов должны совпадать. Для этого

введем функции округления: x – функция 1 2 3 4 5

пола, которая определяется как наибольшее x

целое, меньшее или равное x, а именно Рисунок 1 - Одномерная линейная кусочно-

x = n ⇔ x – 1 < n ≤ x; x – функция потолка, непрерывная функция
© Подгорный С.А., Меретуков З.А.,
Кошевой Е.П., Косачев В.С., 2013

10
Вестник ВГУИТ, №2, 2013
Как видно из представленного графи- вестные функции и операции над ними входят
ка (рисунок 1) функция u(x) на границах во все соотношения задачи только в первой
интервала (a, b) обращается в ноль. Использу- степени, метод конечных элементов получил
ем (1) для описания пробной функции на про- достаточно полное математическое обоснова-
извольном интервале, задаваемом полушири- ние. В дальнейшем используем линейную за-
ной h = (b – a)/2 одномерного единичного ко- дачу, решение которой метод конечных эле-
нечного элемента относительно узла X i . В этом ментов сводит к решению систем линейных
случае параметры кусочно-линейной непре- алгебраических уравнений. Рассмотрим при-
рывной функции (1) соответственно равны менение метода для описания одномерного
a = X i - h, b = X i + h, c = 0, d = 1, а функция температурного профиля, задаваемого в
приобретает вид: начальный момент времени (τ=0 – начальное
условие). В этом случае пробная функция,
x + h − Xi определяемая уравнением (3), представлена
ϕ i (x ) = {Φ[x − ( X i − h )] − Φ(x − X i )}⋅ +
h линейной комбинацией функций (2) с коэффи-
h − x + Xi циентами u i =u i (0):
+ {Φ ( x − X i ) − Φ[x − ( X i + h )]}⋅ , (2)
h
n
Функции φ i (x) изображаются в виде ло- u ( x ) = ∑ ui ⋅ ϕ i ( x ) , (3)
маных и определяются конечным числом па- i =1

раметров - своими узловыми значениями. На

Для того чтобы u(x i ) = u i во всех узлах
графике (рисунок 2) показана функция такого
X i , функции φ i (x) должны удовлетворять усло-
семейства.
виям φ i (X i ) = 1 и φ i (X j ) = 0 для всех узлов X j
2
Xi−h Xi+h при j≠i. Кроме того, чтобы выполнялись гра-
ничные условия первого рода, следует поло-
жить u 0 = u n+1 = 0. Метод конечных элементов
1 d оперирует в качестве φ i (x) кусочно-
полиномиальными функциями, отличными от
ϕ i( x) нуля в пределах небольшого числа элементов
вблизи узла X i . Именно это делает метод мак-
0 c симально
эффективным. Поскольку u(x) по своему физи-
ческому смыслу должна быть непрерывной
функцией, выберем φ i (x) в виде кусочно-
−1
1 2 3 4 5 линейных функций, отличных от нуля на двух
x элементах (рисунок 2). Каждая такая функция
Рисунок 2 – Кусочно-линейная функция (2) для φ i (x), i = 1, 2, …, n, равна единице в X i и нулю
решения континуальной одномерной задачи мето- во всех остальных узлах. При этом набор
дом конечных элементов функций u(x) будет состоять из непрерывных
функций, линейных в пределах элементов с
Метод конечных элементов заменяет за- изломами в узлах и определяемых своими уз-
дачу отыскания функции на задачу отыскания ловыми значениями u i , i = 1, 2, …, n. На кон-
конечного числа ее приближенных значений в цах интервала X min …X max они обращаются в
отдельных точках-узлах. При этом, если ис- нуль. Каждую из таких функций можно изоб-
ходная задача относительно функции состоит разить в виде ломаной линии. Для определения
из дифференциального уравнения с соответ- параметров u i , используемых в уравнении (3),
ствующими граничными условиями, то задача сформируем систему линейных алгебраиче-
метода конечных элементов относительно ее ских уравнений методом Галеркина, интегри-
значений в узлах представляет собой систему руя произведение пробной функции на семей-
алгебраических уравнений. С уменьшением ство кусочно-линейных непрерывных функций
максимального размера элементов увеличива- по области существования решения:
ется число узлов и неизвестных узловых пара-
метров. Вместе с этим, повышается возмож-
∫ [y ( x) ⋅ ϕ (x )]dx, (4)
 
X max n X max

ность более точно удовлетворить уравнениям ∫ ϕ j ( x ) ⋅ ∑ ui ⋅ ϕ i ( x ) dx = j

X min  i =1 
задачи и тем самым приблизиться к искомому
X min

решению. Для линейных задач, когда неиз-

11
Вестник ВГУИТ, №2, 2013
где j=1,2,…,n; y(x) – функция начального тем- Таким образом, уравнение (4) принимает
пературного профиля. матричную форму:
В полученной системе линейных алгеб-
раических уравнений определенные интегралы 2
3
1
6 0 0 0 u1 1
в левой части образуют квадратную матрицу 1
6
2
3
1
6 0 0 u2 1
ленточного типа трехдиагональной структуры, 0 1
6
2
3
1
6 0 ⋅ u3 = 1
образованной двумя видами интегралов на 0 0 1 2 1 u n−1 1
(1)
главной диагонали матрицы:
6 3 6

0 0 0 1
6
2
3 un 1
2 X max − X min
X max

mi ,i = ∫ ϕ (x ) ⋅ ϕ (x )dx = 3 ⋅
i i
n +1
, (5) Здесь матрица является трехдиагональ-
X min
ной, и решение может быть получено методом
Гаусса, матричным методом (умножением на
где i=1,2,…,n. обратную матрицу) или прогонки. Для систем
И над и под главной диагональю: малой размерности метод решения не суще-
ственен, но с увеличением точности решения
− X min
X max

mk ,k +1 = ml ,l −1 = ∫ ϕ (x )⋅ ϕ (x )dx = 6 ⋅
1 X
l −1
max
, (6) необходимо увеличивать число конечных
n +1
l
X min элементов на области существования решения,
а это приводит к увеличению числа неизвест-
где k=1,2,…,n-1; l=2,3,…,n. ных в уравнении (8). Для определения влияния
Упростим задачу, полагая, что y(x)=1. В числа конечных элементов на интервале X min …
этом случае столбец свободных членов p i X max на точность начальной аппроксимации
представляет собой величину, определяемую решали систему (8) матричным методом в ин-
интегралом: женерной среде MathCAD при различных чис-
лах n, определяя среднее значение u(x) на этом
X max
X max − X min интервале. Относительная ошибка при различ-
pi = ∫ ϕ (x ) dx =
i
n +1
, (7 )
ном числе конечных элементов представлена
ниже.
X min

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

элементов матричным методом

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

12
Вестник ВГУИТ, №2, 2013
чения и применим только для систем линей- 1 1 1
α 0 = − , α1 = − ,, α n −1 = − , (9)
ных алгебраических уравнений ленточного 4 4 + α0 4 + α n−2
типа. При применении метода конечных эле-
ментов ширина полосы ленточной матрицы 3 6 − β0 6 − β n−2
зависит от нумерации узлов. В некоторых β0 = , β1 = ,  , β n −1 = , (10)
2 4 + α0 4 + α n−2
случаях исходная постановка задачи может
оказаться настолько плохой, что даже метод
Эти величины используются для расчета
конечных элементов не может помочь. В та-
весовых коэффициентов u i :
ком случае постановку задачи необходимо
менять. При этом имеет место система алгеб- 6 − β n −1
раических уравнений, в которой малые изме- un = , u n −1 =
4 + α n −1
нения коэффициентов или свободных членов
приводят к значительному изменению реше- = α n −1 ⋅ u n + β n −1 ,, u1 = α 0 ⋅ u 2 + β 0 , (11)
ния. Такие системы уравнений носят название
плохо обусловленных. Рассмотрим примене- Для определения влияния числа конеч-
ние метода прогонки для представленной ных элементов на интервале X min …X max на
выше системы. Метод прогонки является точность начальной аппроксимации решали
двухшаговым. Вначале вычисляем вспомога- систему (9), (10) и (11) в инженерной среде
тельные величины α i , β i : MathCAD при различных числах n, определяя
среднее значение u(x) на этом интервале.

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

элементов методом прогонки

Относительная ошибка при различном быть представлена произведением координат-

числе конечных элементов представлена выше ных и временных функций:
(рисунок 4). Из представленных на рисунке
данных видно, что метод конечных элементов Ψ (x ) =
n

позволяет снизить относительную ошибку

∑ u (τ ) ⋅ ϕ ( x ) ,
i =1
i i (12)
начальной аппроксимации практически до
значения менее одного процента при исполь- с граничными условиями Ψ(0) = Ψ(1) = 0.
зовании 40 и более конечных элементов. Ис- Однако поскольку уравнение содержит вто-
пользуя значения весовых коэффициентов u i рую производную по координате, а уже
как значения неизвестных временных функций первая производная конечного элемента
u(τ), можно перейти к решению краевой зада- терпит разрывы непрерывности в узлах,
чи. В этом случае пробная функция может воспользуемся следующим приемом. Обо-
значим R(x) = u(τ)·φ"(x)-u′(τ)·φ(x) невязку ис-

13
Вестник ВГУИТ, №2, 2013
ходного дифференциального уравнения. Точ- структуры, образованной двумя видами инте-
ное решение дает R(x) = 0. Смягчим выполне- гралов на главной диагонали матрицы:
ние этого условия, потребовав, чтобы оно вы-
полнялось только для n функций, которые сов- X max

падают с пробными - u(τ)·φ(x). Такой прием

mi ,i = ∫ ϕ ′(x ) ⋅ ϕ ′(x )dx = n + 1,
X min
i i (15)
носит название метода Галёркина. Выполним
где i=1,2,…,n.
для невязки интегрирование по частям при
И над и под главной диагональю:
условии φ(x) = φ j (x) и φ j (0) = φ j (1) = 0, тогда
получим систему первого порядка, как по вре- X max
n +1
менной составляющей, так и по координатной: mk ,k +1 = ml ,l −1 =
X min
∫ ϕ ′ (x ) ⋅ ϕ ′ (x )dx = − 2 , (16)
l l −1

X max

∫ [u (τ ) ⋅ ϕ ′′(x ) − u ′ (τ ) ⋅ ϕ (x )] ⋅ ϕ (x )dx =
i i i i j
Эта матрица симметричная, что харак-
X min терно для метода Галёркина. Кроме того,
произведение отлично от нуля только при j = i,
∫ [− u (τ ) ⋅ ϕ ′ (x ) ⋅ ϕ ′ (x ) − u ′ (τ ) ⋅ ϕ (x ) ⋅ ϕ (x )]dx = 0 , (13)
X max

− i i j i i j j = i±1, когда соответствующие два элемента,

которые несут на себе пробные функции, пе-
X min

Теперь уже в задачу входит u', что дает рекрываются. Следовательно, матрица в дан-
систему линейных алгебраических уравнений ном случае оказывается трехдиагональной, так
относительно u i вида: как интегрирование совершается только на
двух соседних элементах. Решение получен-
ной системы алгебраических уравнений позво-
∫ [ϕ ′(x ) ⋅ ϕ ′ (x )]dx =
X max

− u i (τ ) ⋅ i j
ляет найти выражение для производных вре-
X min менных составляющих и получить прибли-
женное решение в форме численного интегри-
∫ [ϕ (x ) ⋅ ϕ (x )]dx,
X max

= u i′ (τ ) ⋅ i j (14) рования, например методом Эйлера. Таким

X min образом, для континуальной задачи метод ко-
нечных элементов осуществляет приближен-
Интегралы правой части (14) представ- ный переход к дискретной задаче на основе
лены выражениями (5) и (6), а интегралы ле- соответствующих кусочно-полиномиальных
вой части в полученной системе линейных ал- функций, отличных от нуля на нескольких со-
гебраических уравнений образуют квадратную седних элементах, содержащих узел X i .
матрицу ленточного типа трехдиагональной Решение (в случае Bi = ∝) может быть
представлено следующей системой дифферен-
циальных уравнений первого порядка:

n + 1 − (n +1 2 ) 0 0 0 u1 2
3
1
6 0 0 0 u1′
− ( 2) n + 1 − ( 2)
n + 1 n + 1 0 0 u2 1
6
2
3
1
6 0 0 u2′
(− 1) ⋅ 0 − (n +12 ) n + 1 − (n +12 ) 0 ⋅ u3 = X max − X min ⋅ 0 16 2 3 16 0 ⋅ u3′ (17)
n +1
0 0 − (n +1 2 ) n + 1 − (n +1 2 ) un −1 0 0 1 6 2 3 1 6 un′ −1
0 0 0 − (n +1 2 ) n + 1 un 0 0 0 1 6 2 3 un′

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

u1 (τ k ) u1 (τ k −1 )  n + 1 − (n+1 2 ) 0  u1 (τ k −1 )
−1
0 0 2
3
1
6 0 0 0
 n+1 
u2 (τ k ) u2 (τ k −1 )  − ( 2 ) n + 1 − ( 2 ) 0
n +1 0 1
6
2
3
1
6 0 0  u2 (τ k −1 )
 − 
u3 (τ k ) = u3 (τ k −1 ) −  0 − ( 2 ) n + 1 − ( 2 ) 0 ⋅ u3 (τ k −1 ) ⋅ ∆τ
X X
+ + ⋅ 0 6 3 16 0 ⋅
n 1 n 1 max min 1 2
(18)
 n + 1 
un−1 (τ k ) un−1 (τ k −1 ) 0 0 − (n+1 2 ) n + 1 − (n+1 2 ) 0 0 16 2 3 16 un−1 (τ k −1 )
 
un (τ k ) un (τ k −1 )  0 0 0 − ( 2) n + 1
n +1 0 0 0 16 2 3  un (τ k −1 )
где Δτ – шаг интегрирования системы (17) по времени.

14
Вестник ВГУИТ, №2, 2013
Таким образом, метод конечных элемен- 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 // Stud-
1 Лыков, А. В. Теория теплопроводности ies in mathematical science. 2012. – V. 4. -№ 2.
[Текст] / А. В. Лыков. - М.: Высшая школа, 3 Akin, J. E. Finite element analysis with
1969. error estimators [Text] / J. E. Akin. - Butterworth-
2 Kosachev, V. S. Using rounding function in Heinemann, 2005. – 512 p.
the problems of finite-element analysis [Text] / V. S. 4 Chen, Z. X. Finite element methods and
Kosachev, E. P. Koshevoy, S. A. Podgorny // Stud- their applications [Text] / Z. X. Chen. – Springer,
ies in mathematical science. 2012. – V. 4. -№ 2. 2005. – 424 p.
3 Akin, J. E. Finite element analysis with 5 Solin, P. Partial differential equations and
error estimators [Text] / J. E. Akin. - Butterworth- the finite element method [Text] / P. Solin. –
Heinemann, 2005. – 512 p. Wiley-Interscience, 2006. – 504 p.
4 Chen, Z. X. Finite element methods and 6 Thomee, V. Galerkin finite element
their applications [Text] / Z. X. Chen. – Springer, methods for parabolic problems [Text] /
2005. – 424 p. V. Thomee. – Springer, 2006. – 459 p.

15