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

Московский государственный технический университет 
имени Н.Э. Баумана 
 
 
 
 
А.В. Котович, И.В. Станкевич 

 
РЕШЕНИЕ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ 

МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ 
 
 
Методические указания к решению задач 
 по курсу «Сеточные методы» 

 
 

Под редакцией В.С. Зарубина 

Москва 

Издательство МГТУ им. Н.Э. Баумана 
2010 

УДК 519.3 
ББК 22.161.8 
К73 

Рецензент И.К. Волков 

Котович А.В. 
К73 
Решение задач теплопроводности методом конечных эле- 
ментов : метод. указания к решению задач по курсу «Сеточные 
методы» / А.В. Котович, И.В. Станкевич; под ред. В.С. Зарубина. – М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. – 84, [4] с. : 
ил.  

 

Приведены формулировки стационарных и нестационарных 
задач теплопроводности. Рассмотрены основные особенности построения численного решения этих задач в рамках конечноэлементной технологии. 
Для студентов 4-го курса факультета ФН МГТУ им. Н.Э. Баумана, изучающих курс «Сеточные методы» и выполняющих соответствующую курсовую работу. Могут быть полезны студентам 
старших курсов других факультетов, изучающим численные методы решения краевых задач. 
Рекомендовано Учебно-методической комиссией НУК ФН. 

                                                                                                     УДК 519.3 
ББК 22.161.8 
 
 
Работа выполнена в рамках гранта поддержки ведущих 
научных школ № НШ-4046.2010.8  

 
 МГТУ им. Н.Э. Баумана, 2010 

ВВЕДЕНИЕ 

На первом этапе при решении задач вычислительной термомеханики определяют температурные поля в исследуемых областях. 
Функции, характеризующие температурные поля, являются решениями задач теплопроводности. В настоящее время особое значение придается решению нелинейных стационарных и нестационарных задач теплопроводности. 
В областях, имеющих сложную геометрическую форму, и при 
сравнительно невысоких требованиях к гладкости функций, входящих в формулировку стационарных и нестационарных задач, 
наиболее перспективны численные методы, среди которых продолжительное время лидирующее положение занимает метод конечных элементов (МКЭ) [1–6].  
Широкому и успешному применению МКЭ способствовали его 
основные свойства, например, такие, как естественность, простота, 
доступность, универсальность и высокая технологичность. Метод 
позволяет проводить численный анализ в областях сложной геометрической формы, учитывать особенности граничных условий и 
теплофизических свойств материалов. Характерной особенностью 
МКЭ является прозрачность основных вычислительных процедур, 
дающих возможность эффективно контролировать обработку данных. Кроме того, МКЭ алгоритмически и программно весьма удобен для объединения с современными методами и средствами 
компьютерной графики. 
При решении стационарных задач теплопроводности с помощью МКЭ используют, как правило, интегральную, в частности 
вариационную, формулировку. Существенным преимуществом 
вариационных формулировок является то, что они позволяют не 
только найти решение, но и оценить его погрешность. В этом 
смысле весьма эффективным оказывается применение двойственных вариационных формулировок. Построение и использование 
двойственных вариационных формулировок для получения апо

стериорных оценок точности температурных полей подробно рассмотрены в работах [7–9]. 
Без ограничения общности не удается дать эквивалентную вариационную постановку нестационарным задачам теплопроводности, поэтому технически простым является использование метода 
Роте [10]. В этом случае аппроксимация строится по временному 
переменному, что позволяет перейти от задачи с параболическим 
оператором к последовательности задач эллиптического типа на 
фиксированных временных шагах, каждая из которых часто допускает вариационную постановку. Полученное на предыдущем 
шаге решение рассматривается как начальное на текущем. Однако 
более широкие возможности для решения нестационарных задач 
дает применение метода взвешенных невязок в форме Галеркина 
[5, 7]. Численное решение нестационарных задач с использованием МКЭ, как правило, осуществляется в соответствии с методом 
Галеркина. Решение реализуется в два этапа. На первом этапе с 
помощью процедур МКЭ выполняется дискретизация по пространству, а на втором применяется какая-либо конечноразностная схема на временном отрезке, приводящая к пошаговой 
процедуре интегрирования по времени. 
Если рассматривается нелинейная нестационарная задача без 
линеаризующих процедур, то на каждом шаге по времени придется решать систему нелинейных алгебраических уравнений с помощью итерационных методов. Для того чтобы этого избежать, 
применяют схемы типа предиктор – корректор. Эти схемы на каждом временном шаге требуют решения двух систем линейных алгебраических уравнений. Существенными недостатками использования подобных схем являются общее усложнение алгоритма 
решения и дополнительные затраты оперативной памяти. Эти 
трудности можно обойти, если нестационарную задачу в каждый 
момент времени решать методом простых итераций с явным заданием скорректированных значений коэффициентов уравнения теплопроводности и граничных условий, применяя метод Галеркина 
для построения матричных соотношений МКЭ. Таким образом, в 
каждой точке временного отрезка решение нелинейной нестационарной задачи заменяется последовательностью решений подобных линейных стационарных задач, отличающихся численными 
значениями коэффициентов уравнения теплопроводности и граничных условий. При этом перед очередной итерацией определя

ют численные значения коэффициентов по полученному на предыдущей итерации решению. 
Пошаговое решение нелинейных уравнений эллиптического и 
параболического типов – известный и достаточно широко используемый метод. Наиболее полные результаты получены для эллиптических уравнений. Особую проблему при итерационном решении составляет сходимость. В силу ограниченности объема этот 
вопрос здесь не рассматривается, а заинтересованные читатели 
могут ознакомиться с результатами работ [11, 12]. 
В методических указаниях рассмотрены формулировки нелинейных стационарных и нестационарных задач теплопроводности 
и основные особенности построения численного решения этих задач в рамках конечно-элементной технологии. Особое внимание 
уделено вопросам применения квадратичных изопараметрических 
конечных элементов. Элементы этого типа хорошо аппроксимируют криволинейные границы геометрически сложных областей и 
имеют высокие интерполяционные характеристики. 

1. ПОСТАНОВКА НЕЛИНЕЙНОЙ СТАЦИОНАРНОЙ 
ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ 

Дано трехмерное евклидово пространство 
3
  с произвольно 
выбранной 
прямоугольной 
декартовой 
системой 
координат 



1
2
3
,
,
x x
x
, в которой положение точки фиксировано радиус
вектором 
i i
x

x
e , где 
ie , =1, 3
i
, – единичные орты координатных 

осей; 
ix , =1, 3
i
, – координаты вектора 
.x  Здесь и ниже по повторяющимся латинским индексам проводится суммирование от 1 до 3, 
а по греческим – нет. 
В ограниченной области 
3
G    с кусочно-гладкой границей 
G

 рассматривается нелинейное уравнение теплопроводности 

 




( , ) ,
,
,
0
ij
j
i
V
T T
q
T



x
x
, 
G

x
, 
(1.1) 

с граничными условиями 

 

1
( )
( )
w
S
T
T

x
x , 
1
S
G

 
x
; 
(1.2) 

2
( , ) ,
( , )
i
ij
j
w
S
n
T T
q
T
 

x
x
, 
2
S
G

 
x
; 
(1.3) 

 


3
( , ) ,
( , )
( )
( )
i
ij
j
w
f
S
n
T T
T
T
T
 
 

x
x
x
x
, 
3
S
G

 
x
, (1.4) 

где 
1
2
3
S
S
S
G



, 




1
2
1
3
mes
mes
S
S
S
S






2
3
mes
0.
S
S


 

Здесь 
( )
T x  – температура и 
1
2
3
( ,
, ,
, ,
, )
T
T T
T
T

; 
( , )
ij
T

x
 – ком
поненты тензора теплопроводности, причем 
( , )
( , )
ij
ji
T
T

 
x
x
 

,
1, 3
i j


 
и 
( , )
0
1, 3
T



 
x
 
G
 
x
 
и 


,
,
T
a b
 
 

,
const
a b 
; 
1
2
( , )
( , )
( )
V
q
T
W
T
W
 

x
x
x  – мощность внутренних 
тепловых источников (стоков); 
1( , )
W
T
x
 – мощность внутренних 
тепловых источников (стоков), зависящая от координат вектора x , 
температуры 
( )
T x  и, возможно, ее производных, 
( ),i
T x
 
1, 3
i 
; 

2( )
W
x  – мощность внутренних тепловых источников (стоков), зависящая только от координат вектора x ; 
( )
w
T
x  – заданная температура поверхности 
1
S ; 
( , )
w
q
T
x
 – плотность теплового потока 
в направлении внешней нормали к поверхности 
2
S ; 
( , )
w
T

x
 – 
коэффициент теплоотдачи на поверхности 
3
S ; 
( , )
0
w
T


x
 


3
( , )
,
T
S
a b



x
; 
( )
f
T
x  – температура среды у поверхности 
3
S ; 

in  – компоненты единичного вектора внешней нормали 
i i
n

n
e . 
Запятой с индексом обозначена операция дифференцирования по 
пространственным координатам 

1
2
3
,
,
x x
x
. 
Предположим, что коэффициенты уравнения (1.1) и входящие в 
граничные условия (1.2)–(1.4) являются непрерывными ограниченными функциями. Потребуем, чтобы 
( , )
ij
T

x
, 
1( , )
W
T
x
, 
( , )
w
q
T
x
 и 

( , )
w
T

x
 имели ограниченные производные по 
,
T  а в случае 

1( , )
W
T
x
 – и по 
,i
T  (
1, 3)
i 
, при всех допустимых значениях соответствующих аргументов. 
Решение нелинейной стационарной задачи (1.1)–(1.4) проводят методом последовательных приближений, в соответствии с 
которым его заменяют решением серии линейных задач с коррекцией всех данных, зависящих от температуры. Для линеари

зации задачи (1.1)–(1.4) зафиксируем некоторую температуру 
( )
T  x  
как 
параметр 
в 
коэффициентах 
уравнения 
(1.1) 

( ,
( ))
ij
T 

x
x
, члене 
( ,
( ))
V
q
T 
x
x
, учитывающем мощность внутренних источников (стоков) теплоты, и в граничных условиях 
(1.3)–(1.4) 
( ,
( ))
ij
T 

x
x
, 
( ,
( ))
w
q
T 
x
x
 и 
( ,
( ))
w
T 

x
x
. Линеаризованная задача имеет вид 

 




( ,
) ,
,
,
0
ij
j
i
V
T
T
q
T





x
x
, 
G

x
; 
(1.5) 

 

1
( )
( )
w
S
T
T

x
x , 
1
S
G

 
x
; 
(1.6) 

 

2
( ,
) ,
( ,
)
i
ij
j
w
S
n
T
T
q
T


 

x
x
, 
2
S
G

 
x
; 
(1.7) 

   



3
( ,
) ,
( ,
)
( )
( )
i
ij
j
w
f
S
n
T
T
T
T
T


 
 

x
x
x
x
, 
3
S
G

 
x
; (1.8) 

После решения задачи (1.5)–(1.8) функцию 
( )
T  x  корректиру
ют: 
def
( )
( )
T
T


x
x  и, если сходимость не достигнута, вновь выполняют решение. 
Численное решение каждой линейной задачи (1.5)–(1.8) предполагает построение соответствующего дискретного аналога. Рассмотрим построение дискретного аналога с помощью процедур 
МКЭ, основанных на вариационной формулировке. 

2. ВАРИАЦИОННАЯ ФОРМУЛИРОВКА СТАЦИОНАРНОЙ 
ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ 

Чтобы перейти к вариационной формулировке, необходимо 
построить функционал, стационарная точка которого одновременно является решением линейной задачи (1.5)–(1.8). Рассмотрим 
следующий функционал: 

 
 
Φ
( , )

G
T =
F
T dx

x
, 
(2.1) 

где 
1
2
3
( ,
, ,
, ,
, )
T
T T
T
T

. 

Предположим, что в области G  подынтегральная функция F  
является непрерывной по совокупности ее аргументов, а также 
существуют и непрерывны все ее частные производные до второго 
порядка включительно. Кроме того, потребуем существования и 
непрерывности частных производных первого порядка на границе 
G

 (по крайней мере, на ее части, состоящей из объединения 

2
3
S
S

). 
Необходимым условием существования стационарной точки 
функционала Φ  будет равенство нулю его первой вариации: 

 

 

Φ
,
0
,
i
i
G

F
F
T
T
dx
T
T




















. 
(2.2) 

Применяя в (2.2) вторую формулу Грина к членам, содержа
щим производные 
,
,i

F
T


 и учитывая, что ( , )
(
),
i
i
T
T

 
 и 
= 0
T

 

на 
1
S , получаем 

2
3
Φ
,
0
,
,
i
i
i
i
G
S
S

F
F
F
Tdx
n
Td s
T
T
T





































. 

Отсюда, принимая во внимание произвольность 
,
T

 имеем 

 
,
0;
,
i
i

F
F
T
T












 
(2.3) 

 

2
3
0
,
i
i S
S

F
n
T




. 
(2.4) 

Уравнение (2.3), как известно, является уравнением Эйлера 
(Эйлера – Лагранжа). Теперь, приравнивая по отдельности слагаемые уравнений (1.5) и (2.3), получаем 

 
V
F
q
T

 

, 
(2.5) 



,
, ,
,
i
ij
i
j
i

F
T
T



 





. 
(2.6) 

Аналогично приравнивая по отдельности слагаемые уравнения 
(2.4) и граничные условия (1.7) и (1.8), имеем 

 

2
2
,
,
i
i
ij
j
w
S
i S

F
n
n
T
q
T









, 
(2.7) 

 



3
3
,
,
i
i
ij
j
w
f
S
i S

F
n
n
T
T
T
T





 




. 
(2.8)  

Из соотношений (2.3)–(2.8) следует такая структура функционала (2.1): 

   

2
3

1
1
Φ
,
,
2
.
2
2
ij
i
j
V
w
w
f
G
S
S
T
T T
q T dx
q Tds
T T
T ds




















 (2.9) 

Достаточным условием локального минимума функционала 
(2.9) является его выпуклость в окрестности стационарной точки 
( ).
T  x  В этом случае вторая вариация функционала  

 

3

3
2
2
2

1
Φ
(
, )
(
)
ii
i
w
i
G
S
T
d x
T
d s












 
(2.10) 

должна быть положительной. 
Из физических соображений функции 
( ,
)
T 



x
 (
1, 3,
 
 

суммирования – нет) и 
( ,
)
w
T 

x
 являются строго положительны
ми для любого набора своих аргументов, поэтому 
2Φ
0


. Таким 
образом, вариационную задачу, эквивалентную стационарной задаче (1.5)–(1.8) при фиксированной 
( )
T  x , можно записать так: 

 
Φ( )
min
T 
, T
D

, 
(2.11) 

где 


1

1
2
( )
( ),
.
W
S
D
u
C G
L G
u
T



