close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Московский государственный технический университет
имени Н.Э. Баумана
А.В. Котович, И.В. Станкевич
РЕШЕНИЕ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ
МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
Методические указания к решению задач
по курсу «Сеточные методы»
Под редакцией В.С. Зарубина
Москва
Издательство МГТУ им. Н.Э. Баумана
2010
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
УДК 519.3
ББК 22.161.8
К73
Рецензент И.К. Волков
К73
Котович А.В.
Решение задач теплопроводности методом конечных элементов : метод. указания к решению задач по курсу «Сеточные
методы» / А.В. Котович, И.В. Станкевич; под ред. В.С. Зарубина. – М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. – 84, [4] с. :
ил.
Приведены формулировки стационарных и нестационарных
задач теплопроводности. Рассмотрены основные особенности построения численного решения этих задач в рамках конечноэлементной технологии.
Для студентов 4-го курса факультета ФН МГТУ им. Н.Э. Баумана, изучающих курс «Сеточные методы» и выполняющих соответствующую курсовую работу. Могут быть полезны студентам
старших курсов других факультетов, изучающим численные методы решения краевых задач.
Рекомендовано Учебно-методической комиссией НУК ФН.
УДК 519.3
ББК 22.161.8
Работа выполнена в рамках гранта поддержки ведущих
научных школ № НШ-4046.2010.8
 МГТУ им. Н.Э. Баумана, 2010
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ВВЕДЕНИЕ
На первом этапе при решении задач вычислительной термомеханики определяют температурные поля в исследуемых областях.
Функции, характеризующие температурные поля, являются решениями задач теплопроводности. В настоящее время особое значение придается решению нелинейных стационарных и нестационарных задач теплопроводности.
В областях, имеющих сложную геометрическую форму, и при
сравнительно невысоких требованиях к гладкости функций, входящих в формулировку стационарных и нестационарных задач,
наиболее перспективны численные методы, среди которых продолжительное время лидирующее положение занимает метод конечных элементов (МКЭ) [1–6].
Широкому и успешному применению МКЭ способствовали его
основные свойства, например, такие, как естественность, простота,
доступность, универсальность и высокая технологичность. Метод
позволяет проводить численный анализ в областях сложной геометрической формы, учитывать особенности граничных условий и
теплофизических свойств материалов. Характерной особенностью
МКЭ является прозрачность основных вычислительных процедур,
дающих возможность эффективно контролировать обработку данных. Кроме того, МКЭ алгоритмически и программно весьма удобен для объединения с современными методами и средствами
компьютерной графики.
При решении стационарных задач теплопроводности с помощью МКЭ используют, как правило, интегральную, в частности
вариационную, формулировку. Существенным преимуществом
вариационных формулировок является то, что они позволяют не
только найти решение, но и оценить его погрешность. В этом
смысле весьма эффективным оказывается применение двойственных вариационных формулировок. Построение и использование
двойственных вариационных формулировок для получения апо3
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
стериорных оценок точности температурных полей подробно рассмотрены в работах [7–9].
Без ограничения общности не удается дать эквивалентную вариационную постановку нестационарным задачам теплопроводности, поэтому технически простым является использование метода
Роте [10]. В этом случае аппроксимация строится по временному
переменному, что позволяет перейти от задачи с параболическим
оператором к последовательности задач эллиптического типа на
фиксированных временных шагах, каждая из которых часто допускает вариационную постановку. Полученное на предыдущем
шаге решение рассматривается как начальное на текущем. Однако
более широкие возможности для решения нестационарных задач
дает применение метода взвешенных невязок в форме Галеркина
[5, 7]. Численное решение нестационарных задач с использованием МКЭ, как правило, осуществляется в соответствии с методом
Галеркина. Решение реализуется в два этапа. На первом этапе с
помощью процедур МКЭ выполняется дискретизация по пространству, а на втором применяется какая-либо конечноразностная схема на временном отрезке, приводящая к пошаговой
процедуре интегрирования по времени.
Если рассматривается нелинейная нестационарная задача без
линеаризующих процедур, то на каждом шаге по времени придется решать систему нелинейных алгебраических уравнений с помощью итерационных методов. Для того чтобы этого избежать,
применяют схемы типа предиктор – корректор. Эти схемы на каждом временном шаге требуют решения двух систем линейных алгебраических уравнений. Существенными недостатками использования подобных схем являются общее усложнение алгоритма
решения и дополнительные затраты оперативной памяти. Эти
трудности можно обойти, если нестационарную задачу в каждый
момент времени решать методом простых итераций с явным заданием скорректированных значений коэффициентов уравнения теплопроводности и граничных условий, применяя метод Галеркина
для построения матричных соотношений МКЭ. Таким образом, в
каждой точке временного отрезка решение нелинейной нестационарной задачи заменяется последовательностью решений подобных линейных стационарных задач, отличающихся численными
значениями коэффициентов уравнения теплопроводности и граничных условий. При этом перед очередной итерацией определя4
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ют численные значения коэффициентов по полученному на предыдущей итерации решению.
Пошаговое решение нелинейных уравнений эллиптического и
параболического типов – известный и достаточно широко используемый метод. Наиболее полные результаты получены для эллиптических уравнений. Особую проблему при итерационном решении составляет сходимость. В силу ограниченности объема этот
вопрос здесь не рассматривается, а заинтересованные читатели
могут ознакомиться с результатами работ [11, 12].
В методических указаниях рассмотрены формулировки нелинейных стационарных и нестационарных задач теплопроводности
и основные особенности построения численного решения этих задач в рамках конечно-элементной технологии. Особое внимание
уделено вопросам применения квадратичных изопараметрических
конечных элементов. Элементы этого типа хорошо аппроксимируют криволинейные границы геометрически сложных областей и
имеют высокие интерполяционные характеристики.
1. ПОСТАНОВКА НЕЛИНЕЙНОЙ СТАЦИОНАРНОЙ
ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ
Дано трехмерное евклидово пространство  3 с произвольно
выбранной прямоугольной декартовой системой координат
 x1 , x2 , x3  , в которой положение точки фиксировано радиусвектором x  xi ei , где ei , i = 1, 3 , – единичные орты координатных
осей; xi , i = 1, 3 , – координаты вектора x. Здесь и ниже по повторяющимся латинским индексам проводится суммирование от 1 до 3,
а по греческим – нет.
В ограниченной области G  3 с кусочно-гладкой границей
 G рассматривается нелинейное уравнение теплопроводности
 ij ( x,T ) T , j  , i  qV  x, T   0 ,
x G ,
(1.1)
с граничными условиями
T ( x ) S  Tw ( x ) , x  S1   G ;
1
(1.2)
5
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ni  ij ( x , T )T , j
ni  ij ( x , T )T , j
S3
S2
 qw ( x , T ) , x  S 2   G ;

(1.3)

  w ( x , T ) T ( x )  T f ( x ) , x  S3   G , (1.4)
где S1  S2  S3   G , mes  S1  S 2   mes  S1  S3   mes  S2  S3   0.
Здесь T ( x ) – температура и T  (T , T ,1 , T ,2 , T ,3 ) ; ij ( x , T ) – компоненты тензора теплопроводности, причем ij ( x , T )   ji ( x , T )
i, j  1, 3
и
  ( x , T )  0   1, 3
x  G
и
T   a, b  ,
a, b  const ; qV ( x , T )  W1 ( x, T )  W2 ( x ) – мощность внутренних
тепловых источников (стоков); W1 ( x , T ) – мощность внутренних
тепловых источников (стоков), зависящая от координат вектора x ,
температуры T ( x ) и, возможно, ее производных, T ( x ),i i  1, 3 ;
W2 ( x ) – мощность внутренних тепловых источников (стоков), зависящая только от координат вектора x ; Tw ( x ) – заданная температура поверхности S1 ; qw ( x , T ) – плотность теплового потока
в направлении внешней нормали к поверхности S 2 ;  w ( x , T ) –
коэффициент теплоотдачи на поверхности S3 ;  w ( x , T )  0
( x , T )  S3   a, b  ; T f ( x ) – температура среды у поверхности S3 ;
ni – компоненты единичного вектора внешней нормали n  ni ei .
Запятой с индексом обозначена операция дифференцирования по
пространственным координатам  x1 , x2 , x3  .
Предположим, что коэффициенты уравнения (1.1) и входящие в
граничные условия (1.2)–(1.4) являются непрерывными ограниченными функциями. Потребуем, чтобы ij ( x , T ) , W1 ( x , T ) , qw ( x , T ) и
 w ( x , T ) имели ограниченные производные по T , а в случае
W1 ( x , T ) – и по T ,i (i  1, 3) , при всех допустимых значениях соответствующих аргументов.
Решение нелинейной стационарной задачи (1.1)–(1.4) проводят методом последовательных приближений, в соответствии с
которым его заменяют решением серии линейных задач с коррекцией всех данных, зависящих от температуры. Для линеари6
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
зации задачи (1.1)–(1.4) зафиксируем некоторую температуру
T  ( x ) как параметр в коэффициентах уравнения (1.1)
 ij ( x , T  ( x )) , члене qV ( x , T  ( x )) , учитывающем мощность внутренних источников (стоков) теплоты, и в граничных условиях
(1.3)–(1.4) ij ( x, T  ( x )) , qw ( x , T  ( x )) и  w ( x, T  ( x )) . Линеаризованная задача имеет вид

ij ( x , T




) T , j , i  qV x , T   0 , x  G ;
T ( x ) S  Tw ( x ) , x  S1   G ;
1
 ni  ij ( x , T  )T , j
 ni  ij ( x , T  )T , j
S3
S2
 qw ( x , T  ) , x  S 2   G ;

(1.5)
(1.6)
(1.7)

  w ( x , T  ) T ( x )  T f ( x ) , x  S3   G ; (1.8)
После решения задачи (1.5)–(1.8) функцию T  ( x ) корректируdef
ют: T  ( x )  T ( x ) и, если сходимость не достигнута, вновь выполняют решение.
Численное решение каждой линейной задачи (1.5)–(1.8) предполагает построение соответствующего дискретного аналога. Рассмотрим построение дискретного аналога с помощью процедур
МКЭ, основанных на вариационной формулировке.
2. ВАРИАЦИОННАЯ ФОРМУЛИРОВКА СТАЦИОНАРНОЙ
ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ
Чтобы перейти к вариационной формулировке, необходимо
построить функционал, стационарная точка которого одновременно является решением линейной задачи (1.5)–(1.8). Рассмотрим
следующий функционал:
Φ T  =  F ( x , T ) dx ,
(2.1)
G
где T  (T , T ,1 , T ,2 , T ,3 ) .
7
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Предположим, что в области G подынтегральная функция F
является непрерывной по совокупности ее аргументов, а также
существуют и непрерывны все ее частные производные до второго
порядка включительно. Кроме того, потребуем существования и
непрерывности частных производных первого порядка на границе
 G (по крайней мере, на ее части, состоящей из объединения
S 2  S3 ).
Необходимым условием существования стационарной точки
функционала Φ будет равенство нулю его первой вариации:
 F

F
Φ   
T 
 T ,i   dx  0 .
 T

  T ,i 

G 
(2.2)
Применяя в (2.2) вторую формулу Грина к членам, содержаF
щим производные
, и учитывая, что (T ,i )  (T ),i и T = 0
T ,i
на S1 , получаем
 F  F
Φ   

 T  T ,i
G
 
 F 
 ,i  Tdx    ni
Td s  0 .
T ,i 
 
S2  S3 
Отсюда, принимая во внимание произвольность T , имеем
F  F

T  T ,i
ni
F
T ,i

 ,i  0;

(2.3)
0.
(2.4)
S 2  S3
Уравнение (2.3), как известно, является уравнением Эйлера
(Эйлера – Лагранжа). Теперь, приравнивая по отдельности слагаемые уравнений (1.5) и (2.3), получаем
F
  qV ,
T
8
(2.5)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 F 

 ,i   ijT ,i , j .
 T , i 


(2.6)
Аналогично приравнивая по отдельности слагаемые уравнения
(2.4) и граничные условия (1.7) и (1.8), имеем
ni
ni
F
T , i
F
T , i
S2
  ni  ijT , j  qw  ,
S2

S3

  ni  ijT , j   w T  T f  .

 S3
(2.7)
(2.8)
Из соотношений (2.3)–(2.8) следует такая структура функционала (2.1):
ΦT  
1
1
ij T ,i T , j  2 qV T  dx   qwTds    w  T T f


2G
2
S
S
2
3

 T ds. (2.9)

Достаточным условием локального минимума функционала
(2.9) является его выпуклость в окрестности стационарной точки
T  ( x ). В этом случае вторая вариация функционала
2 Φ  
3
 ii (T ,i )2 d x    w (T )2 d s
G i 1
(2.10)
S3
должна быть положительной.
Из физических соображений функции   ( x , T  ) (   1, 3,
суммирования – нет) и  w ( x , T  ) являются строго положительны-
ми для любого набора своих аргументов, поэтому 2 Φ  0 . Таким
образом, вариационную задачу, эквивалентную стационарной задаче (1.5)–(1.8) при фиксированной T  ( x ) , можно записать так:
Φ (T )  min , T  D ,

(2.11)

где D  u  C1 (G )  L2 (G ) , u S  TW .
1
9
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Из построения функционала (2.9) следует, что функция T 0 ( x ) ,
на которой достигается минимальное значение, одновременно
удовлетворяет задаче (1.5)–(1.8) (при фиксированной функции
T  ( x ) ).
В соответствии с общей логикой итерационного решения исходной нелинейной задачи (1.1)–(1.4) функция T  ( x ) корректируdef
ется: T  ( x )  T 0 ( x ), если сходимость не достигнута, вновь выполняется минимизации функционала (2.9), соответствующего
линейной задаче (1.5) – (1.8).
В работе [7] подробно рассмотрены вопросы построения нелинейных функционалов для различных формулировок стационарных задач теплопроводности, а также условия, обеспечивающие их
экстремальные свойства. В данном случае было показано построение квадратичного функционала. Таким образом, решение последовательности подобных линеаризованных дифференциальных
задач заменяется решением эквивалентных им вариационных задач (2.11).
3. ПОСТАНОВКА НЕЛИНЕЙНОЙ НЕСТАЦИОНАРНОЙ
ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ
Пусть заданы в  3 ограниченная область G с кусочногладкой границей  G и в 1 отрезок  0, t  . В замкнутом цилиндре
  G   0, t    4 , где G  G  G , рассмотрим следующую нестационарную задачу для нелинейного уравнения теплопроводности:
T
ij ( x ,T ) T , j , i  qV  x , T   c( x ,T ) ( x ,T )
, ( x , )  G   0, t  , (3.1)



T ( x ,0)  T0 ( x ) , x  G ,
(3.2)
T ( x, ) S  Tw ( x , ) , x  S1   G ,   0 ,
(3.3)
1
10
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ni  ij ( x , T )T , j
 ni  ij ( x , T )T , j
S2
S3
 qw ( x , , T ) , x  S 2   G ,   0 ,


  w ( x , , T ) T ( x , )  T f ( x , ) ,
(3.4)
(3.5)
x  S3   G ,   0,
где S1  S 2  S3   G , mes  S1  S2   mes  S1  S3   mes  S 2  S3   0 ;
T  (T , T ,1 , T ,2 , T ,3 ) ; Λ  ij ( x ,T ) ei e j , ij ( x,T )   ji ( x ,T ) i, j 
~
 1, 3 и   ( x , T )  0   1, 3 (суммирования нет) x  G и
T   a, b  ,
a, b  const ;
qV ( x , , T )  W1 ( x , , T )  W2 ( x , ) ;
 w ( x , , T )  0 ( x , , T )  S3   0, t    a, b  ; n  ni ei ;  – время;
T0 ( x ) – начальная температура; T ( x , ) – температура; ij ( x , T ) –
компоненты тензора теплопроводности; ei  e j – базисные диады;
qV ( x , , T ) – мощность внутренних тепловых источников (стоков);
W1 ( x , , T ) – мощность внутренних тепловых источников (стоков),
зависящая от вектора x , времени  , температуры T ( x , ) и, возможно, ее производных T ( x, ),i , i  1, 3 ; W2 ( x , ) – мощность тепловых внутренних источников (стоков), зависящая только от координат x и времени  ; c( x , T ) – удельная теплоемкость среды,
занимающей область G ;  ( x , T ) – плотность среды, занимающей
область G ; Tw ( x , ) – температура поверхности S1 ; qw ( x , , T ) –
плотность теплового потока на поверхности S2 ;  w ( x , , T ) – коэффициент теплоотдачи на поверхности S3 ; T f ( x , ) – температура среды у поверхности S3 ; ni – компоненты единичного вектора
внешней нормали n .
Предположим, что коэффициенты уравнения (3.1), начальная
температура T0 ( x ) и коэффициенты, входящие в граничные условия (3.3)–(3.5) являются непрерывными ограниченными функциями. Кроме того, потребуем, чтобы ij ( x , T ) , W1 ( x , , T ) , c( x , T ) ,
 ( x , T ) , qw ( x , , T ) и  w ( x , , T ) имели ограниченные производ11
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ные по T , а в случае W1 ( x , , T ) – и по T ,i (i  1, 3) при всех допустимых значениях соответствующих аргументов.
Начальное условие (3.2) может быть согласовано с граничными условиями (3.3)–(3.5), т. е. подстановка T0 ( x ) вместо T ( x , ) в
(3.3)–(3.5) при   0 не должна приводить к нарушению выполнения граничных условий. При соответствующей гладкости функций, входящих в формулировку нестационарной задачи (3.1) –
(3.5), согласование начального и граничных условий является необходимым для существования классического решения [13, 14].
Данное согласование не считается обязательным требованием существования единственного обобщенного решения, т. е. решения,
удовлетворяющего некоторому интегральному равенству [13]. В
дальнейшем под решением нестационарной задачи (3.1)–(3.5) понимаем именно обобщенное решение.
4. ПОСТРОЕНИЕ МАТРИЧНЫХ СООТНОШЕНИЙ МКЭ
4.1. Понятие конечного элемента
Прежде всего определим понятие конечного элемента: под конечным элементом будем понимать объект (e) , характеризуемый
занимаемой им в евклидовом пространстве  ограниченной замкнутой областью V ( e ) , координатами узлов x (pe) V ( e) и интерполяционными функциями вида
( e ) ( x ) 
M (e)
 N (pe) ( x) (pe)   N (e)  (e)  ,
(4.1)
p 1
где x V ( e ) ; (e) – идентификационная метка конечного элемента; M ( e ) – число узлов конечного элемента (e) ; p – локальный
номер узла конечного элемента (e) ; ( e) ( x ) – интерполируемая
на V ( e )   функция; (pe) – значения функции ( e) ( x ) в узлах
  – вектор-столбец,
конечного элемента (e) : (pe )  ( e) ( x (pe ) ) ; ( e )
составленный из узловых значений (pe ) ; N (pe ) ( x ) – достаточно
12
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
гладкие линейно независимые функции, называемые в дальнейшем функциями формы конечного элемента (e);  N ( e )  – матрица-строка, составленная из функций формы N (pe ) ( x ) .
В зависимости от задачи в качестве евклидова пространства 
могут рассматриваться пространства 1 ,  2 и 3 .
Из (4.1) следуют очевидные свойства функций формы:
1, q  p
, p, q  1, M ( e) , xq( e ) V ( e) ;
1) N (pe ) ( xq( e ) )  
0,
q
p


2)
M (e)

p1
N (pe ) ( x )  1 ,
x V
(e)
, откуда
M (e)
  N (pe) ( x )  ,i  0,
i  1, 3 .
p 1
Аналитическая конструкция функций формы N (pe ) ( x ) зависит от геометрии замкнутой области V ( e ) , числа узлов M ( e ) и
их расположения, т. е. координат узлов x (pe) . В разд. 5.2 рассмотрены конструкции функций формы квадратичных конечных
элементов.
Ограниченные замкнутые области V ( e ) , на которых определяются конечные элементы, представляют собой, как правило, простейшие геометрические фигуры: отрезки линий, треугольники,
четырехугольники, тетраэдры, пирамиды, призмы, параллелепипеды с плоскими или криволинейными гранями (ребрами). В литературе очень часто замкнутые области V ( e ) называют конечными
элементами, при этом учитывают и расположение узлов в областях
V ( e ) , и тип заданных на них интерполяционных функций. В дальнейшем мы также будем отождествлять замкнутую область V ( e ) и
определенный на ней конечный элемент (e) , если это не будет
приводить к противоречию. На рис. 4.1 в качестве примера представлены квадратичные конечные элементы, которые могут быть
использованы для решения задач теплопроводности в геометрически сложных областях разной размерности.
13
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а
б
в
г
е
д
ж
Рис. 4.1. Квадратичные конечные элементы:
а – 3-узловой; б – 6-узловой; в – 8-узловой; г – 10-узловой; д – 13-узловой;
е – 15-узловой; ж – 20-узловой
Аппроксимируем ограниченную замкнутую область G  3 ,
где решаются задачи (1.1)–(1.4) и (3.1)–(3.5), ограниченной
замкнутой областью Gh   3 , состоящей из объединения ко14
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
нечных элементов V ( e ) . При этом должны выполняться следующие условия:
kel
1) Gh =  V ( e ) , где kel – число конечных элементов;

e=1
(i )
2) mes V


 = 0 , i, j  1, k
 G  0 , i  1, kel , V (i )  intV (i ) ;
mes V (i )  V ( j )
el ,
i j;
3) V (i )  V ( j )    S (ij ) , где  – пустое множество; S (ij ) – общая грань (вершина, ребро) для конечных элементов V (i ) ; V ( j ) –
замкнутая область размерности r , здесь r  0 , если S (ij ) – общая
вершина; r  1 , если S (ij ) – общее ребро; r  2 , если S (ij ) – общая
грань;
4) узлы соседних конечных элементов совместны, т. е. если для
элементов с номерами (i ) и ( j ) S (ij )   , то  p, q : x (pi )  xq( j ) , где
p {1, …, M (i ) } и q {1, …, M ( j ) } – локальные номера узлов элементов с номерами (i ) и ( j ) ;
5) узловые параметры соседних конечных элементов совместны, т. е., если для элементов с номерами (i ) и ( j ) p, q :
x (pi )  xq( j ) S (ij ) , где p{1,…, M (i ) } и q{1,…, M ( j ) } , то (pi ) 


(i )
( x (pi ) )

( j)
( xq( j ) )
 (q j )
и x  S
(ij )
 ( x) 
M (i )
 N n(i ) ( x ) (ni) 
n 1
M ( j)
 N m( j ) ( x ) (mj ) .
m 1
Выполнение условия (5) обеспечивает непрерывность интерполяционных соотношений (4.1) при переходе через общую грань
соседних конечных элементов.
Локальная нумерация узлов предполагает нумерацию только в
пределах множества узлов, принадлежащих данному рассматриваемому элементу (e) , в отличие от глобальной нумерации, при
которой нумеруются все узлы конечно-элементной модели, причем совместные узлы элементов нумеруются только один раз. На
рис. 4.1 показана локальная нумерация узлов различных типов конечных элементов. На рис. 4.2 приведен пример глобальной нумерации узлов конечно-элементной модели.
15
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 4.2. Вариант глобальной нумерации узлов
конечно-элементной модели
Если замкнутая область G имеет каноническую форму, то геометрическую аппроксимацию можно построить таким образом,
что G  Gh . В тех случаях, когда G имеет криволинейную границу, из геометрических соображений следует очевидное асимптотическое соотношение


mes  G  Gh  \  G  Gh   0 ,
(4.2)
которое означает, что если max (mesV ( e) )  0 при kel   , то
e 1, kel
равномерное измельчение сетки конечных элементов позволит получить сколь угодно хорошую геометрическую аппроксимацию G
с помощью Gh . Вообще говоря, соотношение (4.2) является избыточным, так как измельчение сетки достаточно проводить только
вблизи криволинейных участков границы  G и совсем необязательно это делать во всей области G . Необходимо отметить, что
во многих работах исследуется влияние изменения области анализа, возникающее вследствие построения сеток. Как показано, например, в работах [15, 16], для квадратичных изопараметрических
элементов с кусочно-полиномиальной аппроксимацией границы
 G произвольной области G , ошибка в аппроксимации точного
16
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
решения u приближенным uh в энергетической норме не превышает O(h 2 ), где h – некоторый средний шаг по пространству.
Граница области  G аппроксимируется внешними гранями
конечных элементов, находящихся непосредственно вблизи аппроксимируемой границы. На рис. 4.3 представлен пример простейшей аппроксимации криволинейной границы. Внешние грани
Рис. 4.3. Вариант аппроксимации границы области
конечных элементов можно рассматривать как отдельные поверхностные конечные элементы [17, 18]. Учитывая, что в выражение
функционала (2.9) входят поверхностные интегралы, такой подход
с вычислительной точки зрения представляется вполне оправданным. На рис. 4.1, а–в представлен общий вид конечных элементов,
которые можно рассматривать как поверхностные. Однако при
исследовании стержневых и оболочечных конструкций указанные
элементы допустимо использовать одновременно в двух качествах: и как поверхностные, и как объемные. В дальнейшем поверхностные конечные элементы будем обозначать S ( e) . Объединение
всех поверхностных конечных элементов образует границу  Gh
k gr
области Gh :  Gh   S ( e) , где k gr – число поверхностных элеe 1
ментов в сетке конечно-элементной модели. Кроме того,
 Gh  S h1  S h 2  S h3
17
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
( S h1 – аппроксимирует поверхность S1 , на которой заданы граничные условия 1-го рода (1.2) или (3.3); S h 2 – аппроксимирует
поверхность S 2 , на которой заданы граничные условия 2-го рода
(1.3) или (3.4); S h3 – аппроксимирует поверхность S3 , на которой
заданы граничные условия 3-го рода (1.4) или (3.5)).
Предполагается, что любой поверхностный конечный элемент S ( e) целиком лежит на какой-либо одной поверхности, за
исключением, может быть, его собственной границы:
int S ( e )  Sh1  S h 2  S h3 . Это же предположение относится и к
объемным конечным элементам. В случае рассмотрения составных
областей аналогично вводят поверхностные конечные элементы на
внутренних поверхностях, на которых заданы граничные условия
4-го рода, учитывающие контактный теплообмен [7].
4.2. Стационарная задача теплопроводности
После построения сетки конечно-элементной модели из (2.9)
имеем
Φh 
1
1
 ij T ,i T , j  2 qV T  d x   qwTd s    w  T T f

2G
2
Sh 2
Sh 3
h

 Td s. (4.3)

Здесь все данные, входящие в формулировку задачи (1.5)–(1.8),
с учетом результатов работ [12, 13] отнесены к области Gh и ее
границе  Gh .
Функционал Φ h в силу свойства аддитивности может быть
представлен в виде
kel
k gr
e 1
e 1
e)
Φ h   Φ (hV
  Φ(hSe ) ,
(4.4)
где
V
18
1  (e) (e) (e)
ij T ,i T , j  2 qV( e) T ( e )  d x, V ( e )  Gh , (4.5)

(e) 2

e)
Φ(hV

Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Φ(hSe )
 qw( e ) T ( e) ds, S ( e)  S h 2 ,
 ( e )
S

1
   (we)  T ( e )  T f( e)  T ( e) ds, S ( e)  S h3 .
 (e)
2

S
(4.6)
С помощь интерполяционного соотношения (4.1) функции T ( e )
и T ( e ) ,i ( i = 1, 3 ), показанные в (4.5) и (4.6), можно выразить через
глобальный вектор узловых температур T  , компонентами которого являются температуры всех узлов сетки конечно-элементной
модели. При этом будет существовать некоторое отличие в представлении T ( e ) для (4.5) и (4.6), связанное с различием в конструкции функций формы N (pe ) ( x ) ( p  1, M ( e) ) для объемных V ( e ) и
поверхностных S ( e) конечных элементов. В общем случае поверхностные конечные элементы имеют меньше функций формы
(так как у них меньше узлов, рис. 4.1, а–в) и конструкция этих
функций существенно проще, поскольку на одну переменную
меньше, что и оправдывает использование поверхностных конечных элементов. Исключением являются ситуации, когда рассматривают одномерные задачи с теплообменом на боковой поверхности и двухмерные задачи с теплообменом на поверхностях в плане
(рис. 4.4).
Таким образом, для объемного интегрирования имеем
T ( e )   NV( e )   aV( e )  T  ,
(4.7)
T ,1( e )   NV( e) ,1 
 ( e )   ( e )  ( e )
(e)
(e)
T ,2    NV ,2   aV  T    B   aV  T  ,
 (e)   (e) 
T ,3   NV ,3 
(4.8)
где  B ( e)  – матрица градиентов объемного конечного элемента
(e) , компонентами которой являются производные функций формы  NV(e ) ,i  ( i  1, 3 );  aV( e )  – матрица геометрических связей
объемного конечного элемента (e) .
19
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а
б
в
Рис. 4.4. Варианты задания граничных условий теплообмена:
а – 3-узловой конечный элемент; б – 6-узловой конечный элемент; в – 6-узловой
конечный элемент
20
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Аналогично для поверхностного интегрирования
T ( e )   N S( e )   aS( e )  T  ,
(4.9)
где  aS( e )  – матрица геометрических связей поверхностного конечного элемента (e) .
Матрица  aS( e )  имеет следующую структуру: число строк
равно локальному числу узлов M S( e ) рассматриваемого конечного
элемента S ( e) , а число столбцов – общему числу узлов сетки конечно-элементной модели kuzl . Строки состоят из нулей, кроме
позиций (в каждой строке по одной), соответствующих глобальным номерам узлов конечного элемента S ( e ), в этих позициях
стоят единицы. Таким образом, в каждой строке – все нули и одна единица. Матрица  aV( e )  имеет аналогичную конструкцию,
отличие заключается только в числе строк – обычно у объемных
конечных элементов их больше, так как больше локальное число
узлов M V( e ) .
(e)
подставить соотношеТеперь если в выражение (4.5) для ΦhV
ния (4.7) и (4.8), получим
1  T  ( e ) T  ( e )  T  ( e )   ( e ) 
(e)  (e)  
T   aV   B   D   B   2 qV  NV   
2
(e)

e)

Φ(hV
V
(4.10)
  aV( e )  T  dx,
где  D ( e)  – матрица, составленная из ij( e) – компонент тензора
теплопроводности Λ ( e )   ij( e) ( x , T  ) ei  e j . Если оси системы ко~
ординат  x1 , x2 , x3  являются главными осями тензора теплопроводности, то матрица имеет вид
21
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
(e)
11

 D (e)    0



 0
0
 (22e)
0
0 

0 .
(e) 
33

Аналогично подстановка в выражение (4.6) для  (hSe ) соотношения (4.9) дает
 (hSe )
 qw( e )  N ( e)   a ( e)  T  ds , S ( e )  Sh 2 ;
2
 S  S 
 ( e )
S
 2



(e)
(e)
(e) 
(e)  1
   w   N S   aS  T   T f  
2

 S3( e )

  N S( e)   aS( e )  T  ds , S3( e)  S h3 .


 

(4.11)

Минимизация функционала (2.9) приводит к уравнению
 h   T 
T
  kel ( e )
 
 T   e1 hV

k gr

e 1

  (hSe)   0 ,
из которого в силу произвольности вариации T  получаем
  kel ( e )
 
 T   e1 hV

k gr

e 1

  (hSe)   0.
(4.12)
В результате подстановки (4.10) и (4.11) в (4.12) и группировки
слагаемых, имеем




( e ) T
( e ) T  ( e )   ( e ) 
(e) 






a
B
D
B
dx
a
T


   V            V   
e 1
 V (e)



kel
22
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»


k gr
T
T
   aS( e)     (we )  N S( e )   N S( e)  ds   aS( e)  T  
 ( e )

e 1
 S3

kel
   aV( e) 
T
e 1
k gr

e 1

V (e)
T
 aS( e) 


T
qV( e )  NV( e )  dx 


( e )  ( e ) T
( e ) ( e )  ( e ) T

q
N
ds    w T f  N S  ds  .
 ( e ) w  S 

S3( e )
 S2

(4.13)
В соответствии с общепринятой терминологией [2–5] введем в
рассмотрение глобальную матрицу теплопроводности  K  и глобальный вектор узловых тепловых усилий  R . Тогда
kel
 K     aV(e) 
T
e 1

k gr

e 1
T
 aS( e ) 




T
(e)
(e)
   B ( e )   D   B dx   aV( e )  
 

 (e) 
 
  
V



( e )  ( e ) T  ( e ) 


N
N
ds   a ( e )  ,
 ( e ) w  S   S    S 
 S3

kel
T
R    aV(e)  
e 1
V (e)
(4.14)
T
qV( e )  NV( e )  dx 


k gr
T
T
T
   aS( e)    qw( e )  N S( e )  ds    (we ) T f( e )  N S( e )  ds . (4.15)
 ( e )

e 1
S3( e )
 S2

Уравнение (4.13), характеризующее тепловое равновесие в узлах сетки конечно-элементной модели (при фиксированной функции T  ( x ) ), примет вид линейного матричного
23
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 K T   R .
(4.16)
Численные методы решения уравнения (4.16) рассмотрены в
работе [19].
4.3. Нестационарная задача теплопроводности
При решении нестационарных задач теплопроводности, как
отмечалось выше, конечно-элементную дискретизацию по пространству целесообразно выполнять на основе метода взвешенных
невязок в форме Галеркина. Для этого поступают следующим образом. Каждому узлу p (в глобальной нумерации) поставим в соответствие финитную функцию N p ( p  1, kuzl , где kuzl – глобальное число узлов сетки конечно-элементной модели). Носителями
функций N p являются объединения областей V ( e ) тех конечных
элементов (e) , которые содержат данный узел p . Кроме того,
предполагается, что функции N p линейно независимы. Тогда температуру и ее производные по координатам и времени на Gh можно интерполировать с помощью следующих соотношений:
T   N T  ;
(4.17)
T , i   N , i T  , i  1, 3 ;
(4.18)
T   N T  ,
(4.19)
где T , T , i , T – интерполированные значения температуры, ее
производных по координатам и времени;  N  – матрица-строка,
составленная из финитных функций N p ( p  1, kuzl ); T  – векторстолбец из узловых значений температуры Tp ( p  1, kuzl ); T  –
вектор-столбец составленный из узловых значений производных
температуры по времени T p ( p  1, kuzl ).
24
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Далее перенесем все члены уравнения (3.1) в левую часть, подставим в нее интерполированные значения температуры и производных и рассмотрим невязку

   


 

R   ij x , T T , j , i  qV x , , T  c x , T  x , T T , (4.20)
где T  (T , T ,1 , T ,2 , T ,3 ) .
Теперь в соответствии с методом Галеркина умножим невязку
R последовательно на функции N p ( p  1, kuzl ), проинтегрируем
по области Gh и результат приравняем к нулю, получим
 RN p dx  
Gh
Gh
   
 ij x , T T , j , i 




 N p ( x ) dx  0,
  q x , , T  c x , T  x , T T 
 V




 

(4.21)
p  1, kuzl .
С учетом соотношений (4.17)–(4.19) выражение (4.21) образует
систему дифференциальных уравнений 1-го порядка относительно
узловых значений температуры (как функций времени). Для ее
решения необходимо добавить начальное условие (3.2), записанное в соответствующем виде.
Перепишем выражение (4.21), опуская для краткости аргументы:
 p dx  0 ,
  ijT , j  , i  qV  N p dx   c TN
Gh
p  1, kuzl . (4.22)
Gh
Вначале рассмотрим в (4.22) второй интеграл, учитывая (4.19),
получим
 p dx  c   N  N p dx T  
 c TN


Gh

Gh
 c   N1 N p N 2 N p  N k
uzl
Gh
N p  dx T  ;
p  1, kuzl .
(4.23)
25
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В соответствии с общепринятой терминологией [2–5] введем в
рассмотрение глобальную матрицу теплоемкости C  , которая на
основании (4.23) может быть записана в виде
C    c  N T  N  dx .
(4.24)
Gh
Выражение (4.24) для вычислений является крайне неудобным.
Интегрирование удобнее проводить по объемам конечных элементов, а результат суммировать. Для этого вместо глобальной матрицы-строки  N  используем локальные матрицы-строки  NV( e)  ,
это позволит выполнить формирование глобальной матрицы теплоемкости по формуле
kel
C     aV(e) 
T
e 1


T
  c ( e ) (e )  NV( e )   NV( e )  dx   aV( e )  . (4.25)

 
 

 (e)
V

Аналогично рассматриваем первый интеграл в (4.22), при этом
при преобразовании дивергентной части используем вторую формулу Грина. В результате получим соотношения, полностью совпадающие с (4.14) и (4.15). Это подтверждает вывод о том, что
процедуры методов Ритца и Галеркина при решении задач с самосопряженными и положительно определенными операторами приводят к одинаковым по конструкции системам алгебраических
уравнений [20, 21].
Таким образом, нестационарная задача (3.1)–(3.5) после дискретизации по пространству (при фиксированной функции T  ( x ) ) сводится к решению задачи Коши для линейного матричного дифференциального уравнения 1-го порядка
C T    K T   R
(4.26)
T  0  T0  ,
(4.27)
с начальным условием
26
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где T0  – проекция функции T0 ( x ) (начального условия (3.2)) на
узлы сетки конечно-элементной модели.
Особенности построения численного решения задачи (4.26)–
(4.27) рассматриваются в разд. 7.
Замечание. Матрицы  aS( e )  и  aV( e )  , присутствующие в выражениях (4.14), (4.15) и (4.25), непосредственно в реальном вычислительном процессе не используют. Соответствующие глобальные
матрицы формируются с помощью специальных алгоритмов, основанных на применении аппарата матриц-указателей (списков
связности). Построение и применение матриц-указателей подробно рассмотрены в работе [22].
5. ИЗОПАРАМЕТРИЧЕСКИЕ ОТОБРАЖЕНИЯ
И ФУНКЦИИ ФОРМЫ КОНЕЧНЫХ ЭЛЕМЕНТОВ
5.1. Построение изопараметрических отображений
Пусть G   3 , поэтому размерность основного пространства –
пространства глобальных координат  x1 , x2 , x3  – везде равна
трем. Конечные элементы, объединение которых образует сетку
конечно-элементной модели, в общем случае являются криволинейными в системе глобальных координат  x1 , x2 , x3  . Интегралы,
входящие в выражения (4.14), (4.15) и (4.25), удобнее вычислять
в локальных декартовых координатах i , i  1, r , где r – размерность пространства локальных координат. Конечные элементы каждого типа в локальных координатах имеют вид строго
фиксированных замкнутых элементарных областей, удобных для
численного интегрирования с помощью квадратурных формул.
Размерность пространства локальных координат зависит от типа
рассматриваемого конечного элемента. Например, для задания координат узлов конечного элемента, представленного на рис. 4.1, а,
в локальной системе координат достаточно одной оси 1 , а координаты узлов конечных элементов, показанных на рис. 4.1, б и в,
в локальной системе задаются двумя координатами  1 ,  2  . При
этом любой конечный элемент в пространстве локальных коор27
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
динат всегда располагается в единичном кубе: i  [ 1, +1] ,
i  1, r , где r {1, 2, 3} .
Положим для определенности r  3 . Переход от локальных координат  1 ,  2 , 3  к глобальным координатам  x1 , x2 , x3  (рис. 5.1)
осуществляется с помощью изопараметрического отображения
[2 – 5].
Рис. 5.1. Изопараметрическое преобразование координат узлов
конечного элемента
Изопараметрическое отображение определяется заданием одних и тех же интерполяционных функций как для преобразования
координат, так и для представления функций, интерполируемых на
конечных элементах (см. (4.1)):
xi (ξ ) 
M (e)
 N (pe) (ξ ) xip(e)   N (e)   xi(e) , i  1, 3 ,
(5.1)
p 1
где M ( e) – число узлов конечного элемента (e) ; xi (ξ ) – i-я координата вектора x (ξ )  xi (ξ ) ei  3 ; xip( e ) – i-я координата узла с
локальным номером p конечного элемента (e) , p  1, M ( e) ;
x 
(e)
i
– вектор-столбец, составленный из i-х координат всех
M ( e) узлов конечного элемента (e) ; N (pe ) (ξ ) – функции формы
28
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
конечного элемента (e) в локальных координатах
 ξ1 , ξ 2 , ξ 3  ;
 N  – матрица-строка, составленная из функций формы
(e)
N (pe ) () , p  1, M ( e) ; ξ  i эi – радиус-вектор точки в пространстве
 1 , 2 , 3  ;
 1 , 2 , 3  .
локальных координат
эi – векторы базиса в локаль-
ных координатах
Таким образом, из (5.1) имеем
 x1 x2 x3    N (e)   x1(e)
x2( e) x3( e )  ,
(5.2)
  – вектор-столбец, составленный из i-х координат
где xi( e )  xi( e)
всех M ( e) узлов конечного элемента (e) .
Преобразование координат, заданное выражением (5.2), должно быть взаимно однозначным, а значит, обратимым. Отсюда, если
пространство локальных координат является трехмерным, определитель матрицы Якоби данного преобразования не должен принимать нулевые значения внутри Vξ( e ) – образа области V ( e ) конечного
элемента
(e)
в
пространстве
локальных
координат
   
– или, во всяком случае, менять знак. Если пространство локальных координат двухмерное, то первая квадратичная форма,
связанная с поверхностным конечным элементом S ( e) , должна быть
положительно определенной (на образе элемента). Если пространство
локальных координат одномерное, рассматриваемый конечный элемент (e) в пространстве глобальных координат  x1 , x2 , x3  должен
иметь форму дуги, принадлежащей гладкой части кривой, при
этом сама кривая может быть кусочно-гладкой. Для выполнения
указанных требований конечные элементы, аппроксимирующие
область G и ее границу  G , должны быть не слишком искривлены и иметь пропорциональные внутренние размеры. Кроме того,
если на ребрах элементов находятся узлы (по одному на ребро), то
они должны быть расположены вблизи середины ребер (в пространстве глобальных координат  x1 , x2 , x3  или, по крайней мере,
в средней их трети) [15, 16].
29
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
5.2. Функции формы конечных элементов
Конструкция функций формы определяется геометрией конечного элемента, общим числом узлов и их расположением. Выбор
типов конечных элементов для решения той или иной задачи зависит от многих факторов: геометрии области анализа G и ее границы  G , граничных условий, общего объема исходных данных,
степени автоматизации их подготовки и т. п. Однако и при этом
выбор типов конечных элементов в значительной степени остается
интуитивным. Так как применение конечно-элементной технологии предполагает создание пакета прикладных программ, пригодного для решения максимально широкого круга задач, то выбор
типов конечных элементов должен обеспечивать аппроксимационную достаточность и универсальность. При этом остается открытым вопрос об оптимальности (по требуемой точности, времени подготовки данных и непосредственно счета) выбранных типов
конечных элементов для решения конкретной задачи теплопроводности.
Для математического моделирования достаточно широкого
круга трехмерных температурных задач вполне пригодны изопараметрические квадратичные конечные элементы, представленные
на рис. 4.1. Элементы этого типа хорошо аппроксимируют криволинейные границы геометрически сложных областей и имеют высокие интерполяционные характеристики [2–5, 8, 15, 16]. Кроме
того, применение изопараметрической техники позволяет удобно,
экономно и с требуемой точностью работать как с исходными
данными, так и с результатами. Ниже приведена структура функций формы конечных элементов, показанных на рис. 4.1.
В локальной системе координат каждому узлу p 
{1, 2, 3, …, M ( e ) } рассматриваемого конечного элемента (e) ставится в соответствие интерполяционная финитная функция
N (pe ) (ξ ) , носителем которой является образ элемента (e) :
supp N (pe ) (ξ )  Vξ( e)   r ,
где r – размерность пространства локальных координат, r 
{1, 2, 3} ; Vξ( e ) – образ области V ( e ) (или S ( e) ) конечного элемента (e) в пространстве локальных координат.
30
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Эти финитные функции N (pe ) (ξ ) называют функциями формы
конечного элемента (e) . Для их построения в зависимости от типа
рассматриваемого конечного элемента используют полные или
усеченные полиномы от локальных координат. В первом случае
конечные элементы принадлежат лагранжевому семейству, а во
втором – сирендиповому [2]. Функции формы элементов лагранжевого семейства строят с помощью интерполяционных полиномов Лагранжа, а структура функций формы элементов сирендипового семейства, как правило, подбирается. В обоих случаях
функции формы должны обеспечивать выполнение соотношений
(4.1), (4.7)–(4.9), (5.2) и при этом, естественно, обладать свойствами 1 и 2, указанными в разд. 4.1. Следуя рекомендациям работ [2–
5, 15, 16], получим следующие конструкции функций формы.
3-узловой конечный элемент (см. рис. 4.1, а, r  1 ):
N (pe ) (ξ ) = N (pe) (1 ) 
1
1  1 p 1 1 p 1 , p  1, 3 ;
2


N (pe) (ξ ) = N (pe ) (1 )  1  12 , p  2,
где 1 p – фиксированные значения локальных координат узлов,
p  1, 3  11 = 1,0, 12 = 0,0, 13 = 1,0  ;  – текущая координата
точки, 1  [1, 1] .
8-узловой конечный элемент (см. рис. 4.1, в, r  2 ):
N (pe ) (ξ ) 
1
1  1 p 1
4

 1  2 p 2  1 p 1  2 p 2  1 ,
N (pe ) (ξ ) 
1
1  1 p 1
2


 1   2 p 1


где ξ p  1 p , 2 p
нат
узлов,

p  1, 3, 5, 7 ;
 1  2 p 2  
   1 p 2 
2
2

, p  2, 4, 6, 8,
– фиксированные значения локальных коорди-
p  1, 8 ;
ξ   1 ,  2 
–
текущие
координаты,
 1 , 2   [1, 1]  [1, 1] .
31
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
6-узловой конечный элемент (см. рис. 4.1, б, r  2 ).
Данный элемент можно получить из 8-узлового, если в нем совместить узлы 1, 2 и 3. Тогда функции формы 6-узлового элемента
строят на основе функций формы 8-узлового с использованием
специальной коррекции, обеспечивающей непрерывность интерполируемых функций в глобальных координатах при переходе через общую границу соседних конечных элементов [5]. Таким образом, 8-узловой конечный элемент рассматривается в качестве
базового для поверхностных элементов. Следовательно,
3
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ )  N ( e) (ξ ) ;
 (pe) (ξ ) ; N ( e ) (ξ )  N
N1( e ) (ξ )   N
5
2
4
3
p 1
 ( e) (ξ )  2N ( e ) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ )  N ( e) (ξ ) ;
N 4( e ) (ξ )  N
5
7
6
 ( e) (ξ ) ,
N 6( e ) (ξ )  N
8
 (pe ) (ξ ) – функции формы 8-узлового конечного элемента,
где N
1
p  1, 8 ; N ( e ) (ξ )  (1  12 ) (1   22 ) – корректирующий член.
8
Как уже отмечалось, присутствие корректирующего члена
N ( e ) (ξ ) обеспечивает непрерывность интерполяционной функции при переходе через общую границу соседних конечных элементов. Аналогичная коррекция используется ниже и в объемных
элементах.
20-узловой конечный элемент (см. рис. 4.1, ж, r  3 ). Для угловых узлов p  1, 3, 5, 7, 13, 15, 17, 19 имеем
1
N (pe ) (ξ )  (1                        .
8
Для узлов p  2, 4, 6, 8, 9, 10, 11, 12, 14, 16, 18, 20 , находящихся
на серединах ребер,
1
N (pe ) (ξ )  (1              
4
                 ,
32
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где ξ p      – фиксированные значения локальных координат узлов, p  1, 20 ; ξ  (1 ,  2 , 3 ) – координаты текущей точки, (1 , 2 , 3 )  [1, 1]  [1, 1]  [ 1, 1] .
20-узловой конечный элемент рассматривается в качестве базового элемента для построения 10-, 13- и 15-узловых элементов
(см. рис. 4.1, г – е).
15-узловой конечный элемент (см. рис. 4.1, е, r  3 ). Объединяя
узлы {1, 2, 3}, {9, 10} и {13, 14, 15} 20-узлового элемента, получаем треугольную призму, функции формы которой имеют вид
 ( e) (ξ )  N
 ( e) (ξ )  N
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ ) ;
N1( e ) (ξ )  N
1
2
3
2
4
 ( e) (ξ )  N ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e ) (ξ )  2N ( e) (ξ ) ;
N 3( e ) (ξ )  N
5
1
4
6
1
 ( e) (ξ )  N ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ ) ;
N 5( e ) (ξ )  N
7
1
6
8
 ( e) (ξ )  N
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ ) ;
N 7( e ) (ξ )  N
9
10
8
11
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ )  N
 ( e) (ξ )  N
 ( e) (ξ ) ;
N 9( e ) (ξ )  N
12
10
13
14
15
(e)
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ )  N ( e) (ξ ) ;
N11
(ξ )  N
16
12
17
2
(e)
 ( e) (ξ )  2N ( e ) (ξ ) ; N ( e) (ξ )  N
 ( e ) (ξ )  N ( e) (ξ ) ;
N13
(ξ )  N
18
2
14
19
2
(e)
 ( e) (ξ ) ,
N15
(ξ )  N
20
 (pe ) (ξ ) – функции форгде ξ  (1 ,  2 , 3 ) – текущие координаты; N
мы 20-узлового элемента, p  1, 20 .
В данном случае корректирующие члены имеют следующую
конструкцию:
N1(e ) (ξ ) 
1
(1  12 ) (1  22 )(1  3 ) ;
16
33
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
N 2(e ) (ξ ) 
1
(1  12 ) (1   22 )(1 + 3 ) .
16
13-узловой конечный элемент (см. рис. 4.1, д, r  3 ). Данный
элемент получится из 20-узлового элемента, если объединить все
узлы какой-либо одной грани, например {1, 2, 3, 4, 5, 6, 7, 8}. Тогда функции формы записывают следующим образом:
8
 ( e) (ξ ) ;
 (pe) (ξ ) ; N ( e ) (ξ )  N
N1( e ) (ξ )   N
2
9
p 1
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ ) ;
N3( e ) (ξ )  N
10
4
11
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e ) (ξ )  N ( e) (ξ )  N ( e ) (ξ ) ;
N5( e ) (ξ )  N
6
13
1
3
12
 ( e ) (ξ )  2N ( e) (ξ ) ;
N 7( e ) (ξ )  N
14
3
 ( e) (ξ )  N ( e ) (ξ )  N ( e ) (ξ ) ;
N8( e ) (ξ )  N
15
1
3
 ( e ) (ξ )  2N ( e) (ξ ) ;
N9( e ) (ξ )  N
16
2
(e)
 ( e ) (ξ )  N ( e) (ξ )  N ( e ) (ξ ) ;
N10
(ξ )  N
17
2
4
(e)
 ( e ) (ξ )  2N ( e) (ξ ) ;
N11
(ξ )  N
18
4
(e)
 ( e ) (ξ )  N ( e) (ξ )  N ( e ) (ξ ) ;
N12
(ξ )  N
19
1
4
(e)
 ( e) (ξ )  2N ( e ) (ξ ) ,
N13
(ξ )  N
20
1
 (pe ) (ξ ) – функции форгде ξ  (1 ,  2 , 3 ) – текущие координаты; N
мы 20-узлового элемента, p  1, 20 .
Корректирующие члены имеют вид
34
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1
(1  1 ) (1  22 )(1  32 ) ;
16
1
N 2(e ) (ξ )  (1  1 ) (1   22 )(1  32 ) ;
16
1
(e)
N 3 (ξ )  (1  12 ) (1   2 )(1  32 ) ;
16
1
N 4( e) (ξ )  (1  1 ) (1  2 )(1  32 ) .
16
N1(e ) (ξ ) 
10-узловой конечный элемент (см. рис. 4.1, г, r  3 ). Этот элемент получится из 20-узлового элемента, если объединить узлы {1,
2, 3, 4, 5, 6, 7, 8}, {9, 10} и {13, 14, 15}. Тогда для построенного
тетраэдра имеем
8
 ( e) (ξ )  N
 ( e) (ξ ) ;
 (pe) (ξ ) ; N ( e ) (ξ )  N
N1( e ) (ξ )   N
2
9
10
p 1
 ( e) (ξ ) ; N ( e ) (ξ )  N
 ( e) (ξ ) ;
N 3( e ) (ξ )  N
11
4
12
 ( e) (ξ )  N
 ( e) (ξ )  N
 ( e) (ξ )  N ( e ) (ξ )  N ( e ) (ξ ) ;
N 5( e ) (ξ )  N
13
14
15
1
2
 ( e) (ξ )  2N ( e ) (ξ ) ;
N 6( e ) (ξ )  N
16
2
 ( e) (ξ )  N ( e) (ξ )  N ( e ) (ξ )  0,5(1   3 )N ( e) (ξ ) ;
N 7( e ) (ξ )  N
17
2
3
4
 (e ) (ξ )  2N ( e) (ξ )  (1   3 )N (e ) (ξ ) ;
N8( e ) (ξ )  N
18
3
4
 ( e) (ξ )  N ( e) (ξ )  N ( e ) (ξ )  0,5(1   3 )N ( e) (ξ ) ;
N9( e ) (ξ )  N
19
1
3
4
(e)
 ( e) (ξ )  2N ( e ) (ξ ) ,
N10
(ξ )  N
20
1
 (pe ) (ξ ) –
где ξ  (1 , 2 , 3 ) – текущие локальные координаты; N
функции формы 20-узлового элемента, p  1, 20 .
Корректирующие члены записывают в виде соотношений
35
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1
(1  1 ) (1  22 )(1  32 ) ;
16
1
N 2(e ) (ξ )  (1  1 ) (1   22 )(1  32 ) ;
16
1
(e)
N 3 (ξ )  (1  12 ) (1  2 )(1  32 ) ;
16
1
N 4(e ) (ξ )  (1  12 ) (1  22 )(1  3 ) .
16
N1(e ) (ξ ) 
6. ОСОБЕННОСТИ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ
МАТРИЧНЫХ СООТНОШЕНИЙ МКЭ
Интегралы, входящие в выражения (4.14), (4.15) и (4.25), вычисляют, как правило, с помощью квадратурных формул Гаусса. В одномерном случае квадратурная формула Гаусса точна для произвольного полинома степени 2n  1 , где n – число точек интегрирования
[23]. Если n малó, то это может привести к ухудшению сходимости
самой схемы МКЭ и (или) утрате точности решения. Однако если n
слишком велико, то существенны временные затраты на выполнение
процедур численного интегрирования. Таким образом, число точек
интегрирования должно быть достаточным: во-первых, для обеспечения сходимости, а во-вторых, для получения решения с требуемой
точностью. Обычно используют схемы интегрирования 2 ×2 ×2   и
3 
 3 3  , т. е. по две или три гауссовых точки в направлении ло-
кальных осей (1 , 2 , 3 ). Как показывает практика, схема интегри-
рования 3 3 3  дает достаточно точные результаты даже для
сильно искаженных (в глобальных осях) элементов. В табл. 6.1 приведены локальные координаты и весовые коэффициенты для схем с
двумя и тремя гауссовыми точками.
Таблица 6.1
Число гауссовых точек
2
3
36
Номера гауссовых точек
Локальные координаты гауссовых точек
Весовые
коэффициенты
1
2
1
2
3
– 0,577350
+ 0,577350
– 0,774597
0,0
+ 0, 774597
1,0
1,0
0,555555
0,888888
0,555555
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
6.1. Интегрирование по объему
Обозначим

 I c(e)  


V
 I (e)  
 BT DB 
T
(e)
c ( e)  ( e)  NV( e)   NV( e )  dx ;

V
(e)
T
 B ( e )   D ( e)   B ( e )  dx ;

 


I  

(e)
qV
V
(e)
T
qV( e)  NV( e)  dx .
(6.1)
(6.2)
(6.3)
Зафиксируем момент времени   t  и некоторое температурное поле T  ( x ) как параметры рассматриваемой задачи. Тогда,
используя аппарат функций формы, компоненты тензора теплопроводности  ij , удельную теплоемкость c, плотность  и мощность внутренних источников (стоков) qV , входящих в выражения
(4.14), (4.15) и (4.25), можно представить в локальной системе координат следующим образом:
 
(6.4)
c (ξ, T  )   NV( e )  c ( e) ;
 
(6.5)
 (ξ, T  )   NV( e)  ( e) ;
 
(6.6)
 
(6.7)
 ij (ξ, T  )   NV( e )   ij( e) ;
qV (ξ, t  , T  )   NV( e )  qV( e) ,
  – вектор, составленный из значений компонент тензора
теплопроводности, взятых в узлах конечного элемента (e) ; c  –
где ij( e )
(e)
37
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
вектор из значений удельной теплоемкости, взятых в узлах конечного элемента (e) ;  ( e) – вектор из значений плотности, взятых
 
 
в узлах конечного элемента (e) ; qV(e ) – вектор из значений мощности внутренних источников (стоков), взятых в узлах конечного
элемента (e) . Здесь радиус-вектор ξ определяет точку, принадлежащую образу конечного элемента в локальной системе координат, т. е. ξ V( e ) .
Вычисление интегралов (6.1)–(6.3) имеет некоторые особенности, зависящие от типа конечного элемента. Ниже рассмотрено
интегрирование для квадратичных конечных элементов, представленных на рис. 4.1.
3-узловой конечный элемент (см. рис. 4.1, а). Для данного конечного элемента элемент объема d x определяется следующим
образом:
d x  A(1 )
3
  xi ,  
1
i 1
2
d 1 ,
(6.8)
где A(1 ) – площадь поперечного сечения.
Производные координат вычисляют по формуле
xi , 1 
 dN1( e)
dxi
 N (je) , 1 xij( e )  
d 1
 d  1
dN 2( e )
d 1
x ( e ) 
( e )   i1 
dN 3  ( e) 
xi 2  ,
d 1   ( e) 
xi 3 
i  1, 3, (6.9)
здесь xij – глобальные координаты узлов элемента, i , j  1, 3 .
Площадь поперечного сечения A(1 ) можно аппроксимировать
с помощью функций формы
 
A( e) (1 )   NV( e )  A( e) ,
(6.10)
  – вектор, составленный из площадей поперечных сече-
где A( e )
ний, отнесенных к узлам элемента.
38
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Матрица градиентов  B ( e)  имеет вид
(e)
(e)

 (e)   (e)
 B    N1 , 1 N 2 , 1 N 3 , 1  .
(6.11)
Матрица  D ( e)  , состоящая из компонент тензора теплопроводности, в данном случае является скаляром:
(e) 
 D ( e)   11

 
.
(6.12)
Таким образом, с учетом (6.8)–(6.12) имеем
n 
T
 I c(e)    c ( e )( e)  NV( e)   NV( e )  A( e)






i 1 

3

k 1
n 
T
 I ( eT)      B ( e)   D ( e)   B ( e )  A( e )
 


 B DB  i 1  


I    q
( e)
qV
n
i 1 
(e)
V
T
 NV( e )  A( e )


3
  xk , 
k 1
xk ,1
3

2
  xk , 
k 1
2
1
2
1



Hi ;
1 1i

H i ; (6.13)


1 1i

Hi ,


1 1i
где n – число гауссовых точек; i – номер гауссовой точки; 1i –
координата гауссовой точки в локальной системе координат; H i –
весовой коэффициент квадратурной формулы.
6- и 8-узловые конечные элементы (рис. 4.1, б, в). Элемент объема для этого типа конечных элементов определяется по формуле
d x  h(ξ ) EG  F 2 d 1 d  2 ,
(6.14)
где h(ξ ) – толщина конечного элемента.
39
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Коэффициенты квадратичной формы E , G и F задаются соотношениями
3
3
3
i 1
i 1
i 1
E   xi2 ,1 ; G   xi2 , 2 ; F   xi ,1 xi , 2 ,
(6.15)
а частные производные координат вычисляют по формуле, аналогичной (6.9).
Используя функции формы, можно аппроксимировать толщину:
 
h( e) (ξ )   NV( e )  h( e) ,
 
где h( e)
(6.16)
– вектор, составленный из значений толщины конечно-
го элемента, отнесенной к его узлам.
Матрица градиентов  B ( e)  , входящая в выражение для
 I ( eT)  , для данного типа элементов состоит из двух строк:
 B DB 
 N1( e ) ,  ... N ( e)( e ) ,  
1
1
M


 B (e)   
;


 N ( e) , ... N ( e) , 
M ( e ) 2 
 1 2
(6.17)
где M ( e) – число узлов конечного элемента (равно шести или восьми).
Матрица  D ( e)  в данном случае имеет размерность 2  2 :
 ( e)
 D ( e)    11


 0
0 
.
 (22e) 
Окончательно выражения для  I c(e)  ,  I ( eT)  и  I q( e)  при V 
 B DB 
мут вид
40
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
n m
T
I c(e )    c ( e ) ( e ) NV( e )  NV( e )  h( e ) EG  F 2 


 

 


i1 j 1
Hi H j ;
( 1i , 2 j )
n m
T
I ( eT)    B ( e )  D ( e ) B ( e )  h ( e ) EG  F 2 

 


 B DB  i1 j1
H i H j ; (6.18)
( 1i , 2 j )
n m
T
I ( e )    qV( e) NV( e)  NV( e )  h( e ) EG  F 2 
q



 

 V
i1 j 1
Hi H j ,
( 1i , 2 j )
где n , m – число гауссовых точек в направлении локальных осей
(1 ,  2 ) , обычно n  m ; (1i ,  2 j ) – координаты гауссовой точки
(i, j ) в локальной системе координат; H i , H j – весовые коэффициенты квадратурной формулы.
10-, 13-, 15- и 20-узловые конечные элементы. Для указанных
конечных элементов (см. рис. 4.1, г–ж) при замене переменных
элемент объема записывают следующим образом:
d x  d x1d x2 d x3  det J d 1d  2 d 3  det J d  ,
(6.19)
где det J – определитель матрицы Якоби.
Матрица Якоби имеет размерность 3  3 и определяется соотношением
J (e)
 N ( e) ,  ... N ( e )( e ) ,  
1
1 
M
 1


  N1(e ) , 2 ... N ( e )( e ) , 2   x1( e) x2( e ) x3( e )  ,
M




 (e)

(e)
N , ... N ( e ) , 3 
M
 1 3

(6.20)
где M ( e ) – число узлов конечного элемента.
41
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Элементами матрицы градиентов  B ( e )  являются производные функций формы по переменным глобальной системы координат ( x1 , x2 , x3 ) . Учитывая (6.20), имеем
 N ( e ) , 1 ... N ( e )( e ) ,1 
M
 1



(
e
)
(
e
)
(
e
)
 B    N1 , 2 ... N ( e ) ,2   J ( e )


M


 N ( e ) , ... N ( e ) , 
 1 3

M (e) 3 
 
 N ( e) ,  ... N (e )( e ) ,  
1
1 
M
 1

1  ( e )
(e)
 N1 , 2 ... N M ( e ) , 2  . (6.21)


 N ( e) , ... N ( e ) , 
 1 3

M ( e ) 3 
Матрица  D ( e)  для данного типа конечных элементов имеет
размерность 3  3 :
(e)
11

 D (e)    0



 0
0
 (22e)
0
0 

0 .
(e) 
33

Таким образом, с учетом (6.19)–(6.21) для интегралов по объему имеем следующие выражения:
l
m n
 ( e ) ( e )  ( e ) T  ( e ) 
(e) 
 (e) 
Ic  c  NV  NV det J 
i1 j 1 k 1
l
Hi H j H k 
( 1i , 2 j , 3k )
m n
T
I ( eT)  B( e )  D( e ) B( e ) det J ( e ) 

 
 B DB  i1 j1 k 1  
l
( 1i , 2 j , 3k )
m n
T
I ( e )  qV( e ) NV( e )  NV( e ) det J ( e ) 
q
 

 

 V
i1 j 1 k 1
42
Hi H j H k  (6.22)
Hi H j H k 
( 1i , 2 j , 3k )
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где l , m , n – число гауссовых точек в направлении локальных
осей (1 , 2 ,  3 ) , обычно l  m  n ; (1i ,  2 j , 3k ) – координаты
гауссовой точки (i, j , k ) в локальной системе координат; H i , H j ,
H k – весовые коэффициенты квадратурной формулы.
6.2. Интегрирование по поверхности
Обозначим
 I (e)  
 w 
I  
( e)
qw
 (we)  N S( e)   N S( e)  ds;

qw( e )  N S( e )  ds;
S3( e )
T
S2( e )
I  
(e)
 wT f
T


S3( e )
(6.23)
T
 (we) T f  N S( e)  ds.
Для интерполяции (в локальных координатах) на рассматриваемом поверхностном элементе (e) коэффициента теплоотдачи
 (we ) ( x , t  , T  ) , температуры окружающей среды T f( e ) ( x , t  ) или
плотности теплового потока qw( e ) ( x , t  , T  ) используют функции
формы
 
 (we) (ξ, t  , T  )   N S( e)   (we ) ;
(6.24)
 
T f( e) (ξ, t  )   N S( e )  T f( e) ;
(6.25)
 
(6.26)
qw( e ) (ξ, t  , T  )   N S( e)  qw( e) ,
 
где  (we ) – вектор, составленный из значений коэффициента теп-
 
лоотдачи, взятых в узлах конечного элемента (e) ; T f( e ) – вектор,
43
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
составленный из значений температуры окружающей среды, взя-
 
тых в узлах конечного элемента (e) ; qw( e )
– вектор, составлен-
ный из значений плотности теплового потока, взятых в узлах конечного элемента (e) .
Ниже рассмотрены некоторые особенности интегрирования по
поверхности, связанные с геометрией области анализа.
3-узловой поверхностный конечный элемент. При исследовании температурных полей стержневых элементов конструкций
граничные условия по теплообмену могут быть заданы на торцевых и боковой поверхностях (см. рис. 4.1, а и 4.4, а). При задании
граничных условий на торцевых поверхностях интегралы (6.23)
являются скалярными величинами и вычисляются точно:
I  w   (3)
w S3 ;
(3)
I  wT f   (3)
w T f S3 ;
(6.27)
I qw  qw(2) S 2 ,
где S 2 и S3 – площади торцевых поверхностей (см. рис. 4.4, а);
(3)
 (3)
и qw(2) – параметры граничных условий на соответстw , Tf
вующих торцевых поверхностях.
Значение I  w суммируется со значением элемента главной
диагонали матрицы  K  , номер которого равен глобальному номеру узла, размещенному непосредственно на рассматриваемой
торцевой поверхности стержневой системы. Аналогичным образом значения I  wT f и I qw подставляют в вектор  R .
При интегрировании по боковой поверхности элемент площади
d s  P(1 )
3
  xi ,  
i 1
1
2
d 1 ,
(6.28)
где P(1 ) – периметр боковой поверхности конечного элемента (e) .
Используя функции формы, периметр P (1 ) можно интерполировать на конечном элементе следующим образом:
44
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 
P ( e) (1 )   N S( e )  P ( e) ,
(6.29)
  – вектор, составленный из значений периметра боко-
здесь P ( e)
вой поверхности, взятых в узлах конечного элемента (e) .
Таким образом, квадратурные формулы для вычисления интегралов (6.23) по боковой поверхности имеют вид
n
T
 I ( e)     (we)  N S( e)   N S( e)  P ( e)

 

 w 
i 1
3
  xi ,  
k 1
n
 I ( e )    ( e )T ( e )  N ( e )  T  N ( e )  P ( e )
  wT f   w f  S   S 
i 1
n
T
 I ( e)    qw( e )  N S( e )   N S( e )  P ( e )

 

 qw 
i 1
Hi ;
1 1i
3
  xk ,  
3
2
H i ; (6.30)
1
k 1
  xi ,  
k 1
2
1
1 1i
2
Hi ,
1
1 1i
где n – число гауссовых точек; i – номер гауссовой точки; 1i –
локальная координата гауссовой точки; H i – весовой коэффициент квадратурной формулы.
3-узловые поверхностные конечные элементы могут быть использованы при рассмотрении граничных условий на торцевой
поверхности оболочки, если оболочка представляет собой поверхность с краем (рис. 6.1). В этом случае в качестве периметра P(1 )
задается толщина оболочки h(1 ) , отнесенная к узлам на срединной линии торцевой поверхности, а численное интегрирование
граничных условий проводится по формулам (6.30).
6- и 8-узловые поверхностные конечные элементы. Конечные
элементы этого типа (см. рис. 4.1,б, в) используют при моделировании температурных полей оболочек и трехмерных тел. Элемент
поверхности d s записывают следующим образом:
ds 
EG  F 2 d 1 d  2 ,
(6.31)
45
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где коэффициенты 1-й квадратичной формы E , G и F задают
соотношениями (6.15). Отсюда соответствующие квадратурные
Рис. 6.1. Численное интегрирование граничных условий теплообмена
по боковой поверхности оболочки
формулы для вычисления интегралов по поверхности приобретают вид
n
m
T
 I ( e)      (we )  N S( e)   N S( e)  EG  F 2 




 

 w
i 1 j 1

I ( e)T
w f

n
n
 ( 1i , 2 j )
m
T
    (we )T f( e)  N S(e )  EG  F 2 

i 1 j 1 
I    q
(e)
qw
Hi H j ;
m
i 1 j 1
(e)
w
 N S( e ) 


T
EG  F 2 

Hi H j ;
(6.32)
 ( 1i , 2 j )
Hi H j ,
 ( 1i , 2 j )
где n , m – число гауссовых точек в направлении локальных осей
(1 ,  2 ), обычно n  m ; (1i ,  2 j ) – координаты гауссовой точки
46
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
(i, j ) в локальной системе координат; H i , H j – весовые коэффициенты квадратурной формулы.
7. ОСОБЕННОСТИ ЧИСЛЕННОГО РЕШЕНИЯ
ЗАДАЧИ КОШИ
Уравнение (4.26) совместно с начальным условием (4.27) составляют задачу Коши. Для ее решения существуют разные подходы [2–5, 15, 21, 23], однако наибольшее распространение получили два. Первый состоит в том, что производную по времени в
уравнении (4.26) заменяют каким-либо конечно-разностным аналогом, а второй заключается в использовании конечных элементов
во временной области (метод Галеркина).
7.1. Двухслойные разностные схемы
Рассмотрим основные этапы построения разностного аналога
задачи Коши (4.26) и (4.27) в виде семейства двухслойных разностных схем. В пределах временного шага h  n  n 1 векторы
узловых температур T  и правой части  R представим в виде
следующих линейных комбинаций:
T ()  (1  )T n1   T n ;
(7.1)
R()  (1  )Rn1   Rn ,
(7.2)
где  [n 1 , n ]  [0, t ] ,  – весовой множитель   [0, 1] . Здесь
и везде далее векторы T n 1 ,  Rn 1 и T n ,  Rn отнесены к мо-
ментам времени соответственно   n1 и   n .
Разностную аппроксимацию производной по времени можно
представить следующим образом:
 T  T n  T n 1
.

T   
h
(7.3)

Теперь, полагая, что коэффициенты в уравнении (4.26) постоянны на отрезке [n 1 , n ] , и подставляя (7.1), (7.2) и (7.3) в (4.26),
47
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
получим, после очевидных преобразований общее выражение для
двухслойной схемы с весами
C   h  K T n 


 C   h (  1)  K T n 1  h (1  ) Rn 1   Rn . (7.4)
Как известно, фиксированное численное значение параметра
 определяет тип конкретной разностной схемы [15, 21]. Напри-
мер,   0 – схема с разностью вперед,   1 2 – схема Кранка –
Николсона,   2 3 – схема Галеркина,   1 – схема с разностью
назад. Вообще, параметр  может принимать любые численные
значения из отрезке  0, 1 . Схема Кранка – Николсона имеет определенные преимущества перед остальными схемами, поскольку
аппроксимирует уравнение (4.26) по переменной  с порядком
O (h2 ), а остальные указанные схемы имеют более низкий порядок
аппроксимации – O (h ) .
7.2. Трехслойные разностные схемы
Кроме двухслойных схем существуют и используются трехслойные [3, 23]. Рассмотрим на оси времени два фиксированных
отрезка [n 1 , n ] и [n , n1 ] , причем [n 1 , n ]  [n , n 1 ] 
 [n 1 , n 1 ]  [0, t ] и n  n 1  n 1  n  h . Аппроксимируем на
объединенном отрезке [n 1 , n 1 ] производную по времени следующим образом:
 T  T n 1  T n 1

,
T   
2h
(7.5)

затем усредним по узлам вектор температуры:
T  
T n1  T n  T n1
3
.
(7.6)
Здесь и далее вектор T n 1 отнесен к моменту времени   n1 .
Зафиксируем вектор правой части
48
R
в средней точке отрезка
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
[n 1 , n 1 ] и подставим выражения (7.5) и (7.6) в уравнение (4.26),
тогда после преобразований получим
2
2




 C   h  K   T n 1   C   h  K   T n 1 
3
3





2
h  K T n  2h Rn .
3
(7.7)
Другим способом построения разностных схем является метод
Галеркина. Рассмотрим на отрезке [n 1 , n 1 ] три линейные функции
формы, определяющие конечные элементы во временной области:
   n 1
, n 1    n ;
1 
h
N n 1  
0,
 [n 1 , n ] ;

   n 1
n 1    n ;
,
 h



N n  1    n 1 ,      ;
n
n 1
h


0,
 [n 1 , n 1 ] ;

(7.8)
   n
, n    n 1 ;

N n 1   h
 0,
 [n , n 1 ] 

Графики функций формы N n 1 , N n и N n 1 представлены на
рис. 7.1, а. В соответствии с методом Галеркина и учитывая конструкцию функций (7.8), имеем
n 1
 C T    K T   R N n dt  0 ,
(7.9)
n 1
где 0 – нулевой вектор.
49
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
а
б
Рис. 7.1. Графики функций формы для интерполяции по оси времени
Используя функции формы (7.8), векторы производных T  ,
температур T  и правых частей  R можно на отрезке интегрирования [n 1 , n 1 ] представить в виде следующих интерполяционных соотношений:
dN
n
T n  n1 T n1 ;
T   dNd n1 T n1  dN
d
d
T   N n1 T n1  N n T n  N n1 T n1 ;
R  N n1 Rn1  N n Rn  N n1 Rn1 .
50
(7.10)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Здесь и далее вектор  Rn 1 отнесен к моменту времени   n1 .
Если интерполяционные соотношения (7.10) подставить в уравнение (7.9), выполнить интегрирование и необходимые алгебраические преобразования, то получим
1
1
4




 C   h  K   T n 1   C   h  K   T n 1  h  K T n 
3
3
3




1
 h  Rn 1  4 Rn   Rn 1 
(7.11)
3


Аналогично в качестве функций формы N n 1 , N n и N n 1 можно
рассмотреть квадратичные функции, которые на отрезке [ h ,  h ]
задают выражениями
 2  h 
 h2  2
,
[
,
];
,  [h ,  h ];
h
h








N n 1   2h2
N n   h2


 [h ,  h ];
 [h ,  h ];
0,
0,
(7.12)
   h 
,   [ h ,  h ];

N n 1   2h2

 [h ,  h ].
0,
2
Подставляя функции формы (7.12) в (7.10), а затем выведенные
интерполяционные соотношения в уравнение (4.26), после преобразований получим семейство разностных схем
(2  h ) C   (2  h )  K  T  
n 1


  (h  2) C   (h   2 )  K  T n 1 
  4 C   2(2  h2 )  K  T n  (2  h )  Rn 1 
 2 (h2  2 )  Rn  (2  h )  Rn 1 
(7.13)
51
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Фиксируя значение  в диапазоне h     h , из общего выражения (7.13) найдем конкретную схему. Рассмотрим следующие
три варианта:   h ,   0 и    h , тогда из (7.13) соответственно получим
C  T n1  (2 h )  K   3 C  T n1  4 C  T n  2 h Rn1 , (7.14)
C  T n1  C  T n1  2 h  K  T n  2 h Rn ,
 3 C   2h  K T n1  C   4 T n  T n1   2 h Rn1 .
(7.15)
(7.16)
Замечание. Если вместо отрезка [ h , h ] рассмотреть [n1 ,n1 ] ,
то квадратичные функции формы N n 1 , N n и N n 1 будут иметь
более сложный аналитический вид:
 2  (n  n 1 )  n n 1
,  [n 1 , n 1 ];

N n 1  
2h2

 [n 1 , n 1 ];
0,
 2  (n 1  n 1 )  n 1n 1
,  [n 1 , n 1 ];

Nn  
h2

  [n 1 , n 1 ];
0,
(7.17)
 2  (n 1  n )  n 1n
,  [n 1 , n 1 ];

N n 1  
2h2

 [n 1 , n 1 ] 
0,
Графики функций (7.17) представлены на рис. 7.1, б. Соответственно более сложными будут и алгебраические преобразования.
Тем не менее после выполнения всех преобразований из общего
выражения для фиксированных значений   n1 ,   n и   n1
получим те же самые разностные схемы, заданные соотношениями
(7.14)–(7.16).
52
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Еще одну трехслойную схему можно построить, если квадратичные функции формы (7.12) или (7.17) использовать в процедурах метода Галеркина. Подставляя (7.12) в (7.10) и найденные соотношения в (7.9), после интегрирования и необходимых
преобразований окончательно получим
1
1




 C   h  K   T n 1   C   h  K   T n 1 
5
5





8
1
h  K T n  h
5
5
Rn1  8Rn  Rn1  
(7.18)
Рассмотренные выше трехслойные разностные схемы используют начиная со второго шага, а на первом шаге применяется какая-либо двухслойная схема (7.4).
7.3. Диагонализация матрицы теплоемкости
При решении нестационарных температурных задач одной
из острых проблем является численная устойчивость применяемой разностной схемы. Существует ряд определений устойчивости и соответствующих подходов к ее исследованию [15, 21,
23]. Например, двухслойные разностные схемы (7.4) считаются
условно устойчивыми, если   1 2 , и абсолютно устойчивыми,
если   1 2 . Однако устойчивость разностной схемы совсем не
гарантирует отсутствия численных колебаний (осцилляций) узловых значений температуры, особенно заметных на первых
временных шагах при импульсном изменении граничных условий [6, 24].
Для подавления численных колебаний используют следующие
два способа. Первый заключается в том, что, исходя из конкретных условий решаемой температурной задачи, строят аналитическое соотношение между временным h и пространственным hG
шагами [5, 15, 23]. При этом, как правило, предполагается совместное изменение (чаще всего – уменьшение) пространственных
размеров конечных элементов и временного шага h . Ограничение
только временного шага h может оказаться недостаточным для
подавления осцилляций.
53
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Второй способ подавления численных колебаний узловых значений температуры состоит в диагонализации матрицы теплоемкости C  и применении явной схемы интегрирования по временной переменной. Алгоритмам диагонализации посвящено
достаточно много работ [24, 25]. В большинстве обоснование диагонализации носит экспериментальный вычислительный характер.
Треугольные конечные элементы с линейными функциями формы
рассмотрены в работе [25]. Для конечных элементов высокого порядка аппроксимации диагонализацию матрицы теплоемкости C 
можно выполнить, используя для численного интегрирования
(4.25) квадратурные формулы, точки интегрирования которых
расположены непосредственно в узлах конечных элементов [21,
25]. Однако точность их при одинаковом числе точек интегрирования существенно ниже, чем у квадратурных формул Гаусса. Из
приведенных выше разностных схем диагонализация может быть
применена к двухслойной схеме (7.4) при значении параметра
  0 и к трехслойным схемам (7.14) и (7.15).
Диагонализацию глобальной матрицы теплоемкости C  выполняют, используя для ее формирования диагонализированные
локальные матрицы теплоемкости – матрицы теплоемкости конечных элементов C (e )  . Диагонализацию матриц C ( e)  осуществляют по-разному [24, 25]. Укажем два формальных способа. В соответствии с первым диагонализация выполняется по
формулам
 (e) 
C

M (e)
 C(e) ; C (e)  0 , где    и ,   1, M (e) ,
(7.19)
1
(e)
– элементы
здесь M ( e) – число узлов конечного элемента (e) ; C
 ( e) – элементы диагонализированматрицы теплоемкости (4.25); C

ной матрицы теплоемкости. Напомним, что по повторяющимся
греческим индексам суммирования нет.
Для диагонализации по второму способу используют соотношения
54
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 (e) 
C

здесь b 
M (e) M (e)

1 1
(e)
C
и g
b (e)
C ,
g
(7.20)
M (e)
(e)
;
 C
1
 ( e)  0 , где    и ,   1, M (e ) .
C

В первом случае глобальная матрица теплоемкости называется
конденсированной, а во втором – диагонализированной.
8. МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ
АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
8.1. Основные понятия теории итерационных методов
Для решения систем линейных алгебраических уравнений
(СЛАУ) создано много прямых и итерационных методов, некоторые из которых постоянно совершенствуются. Выбор метода решения заданного класса сеточных уравнений, его алгоритмическое
построение и программная реализация в значительной степени зависят от типа решаемой краевой задачи (одномерная или многомерная, стационарная или нестационарная, линейная или нелинейная) и технических характеристик доступной вычислительной
техники. При решении двух- и трехмерных задач большое значение имеет вид области анализа.
В настоящее время известно много работ, посвященных применению прямых методов решения СЛАУ с разреженными матрицами, достаточно полную библиографию по прямым методам
можно найти в работах [26–36]. Прежде всего это относится к развитию семейства метода Холецкого и методов факторизации. Перспективным направлением в разработке прямых методов является
реализация алгоритма быстрого дискретного преобразования Фурье. Наиболее значительные результаты применения прямых методов получены при решении СЛАУ не очень большого порядка в
геометрически простых областях, например типа прямоугольника
и параллелепипеда.
55
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
При решении краевых задач в областях общего вида прямые методы требуют значительных временных затрат и, кроме того, обостряется проблема накопления ошибок округления. В отдельных случаях решение, полученное с помощью прямого метода, подвергается
итерационному уточнению. Достаточно сложна для прямых методов
проблема рационального использования оперативной памяти. Известно, что прямые методы требуют в общем случае больше оперативной памяти по сравнению с итерационными, что связано с особенностями процедуры исключения (проблема восполнения).
Несмотря на успехи в создании и развитии эффективных прямых
методов, в настоящее время решение СЛАУ с разреженными матрицами высокого порядка (десятки и сотни тысяч и даже несколько
миллионов уравнений) и областей сложной формы с помощью итерационных методов представляется более перспективным [37 – 39]. В
пользу выбора итерационных методов говорят такие факторы, как
разреженность матрицы (в процедурах хранятся и используются
только ненулевые элементы нижней или верхней треугольной матрицы), отсутствие затруднений с нумерацией (перенумерацией) узлов
сетки (структура матрицы не обязательно ленточная), простота алгоритмической и программной реализаций, возможность оперативного
контроля достигнутой точности в процессе работы.
Ограничимся рассмотрением некоторых итерационных методов решения сеточных уравнений, которые нашли применение в
конечно-элементных комплексах и пакетах прикладных программ.
Пусть в результате построения дискретного аналога некоторой
краевой задачи с помощью процедур МКЭ получена система n
линейных алгебраических уравнений (см. (4.16), (7.7), (7.11))
Ax  f ,
(8.1)
где x – вектор-столбец узловых неизвестных; f – вектор-столбец
правой части; A  AT  0 . Неравенство означает, что  Ax, x   0
для всех x  0 .
Уравнение (8.1) можно рассматривать как операторное с линейным самосопряженным положительно определенным оператором A , действующим в евклидовом пространстве  n :
A : x  f  im A   n , x   n .
56
(8.2)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для решения уравнения (8.1) применяют двухслойные
B
xk 1  xk
 A xk  f
k 1
(8.3)
и трехслойные
Bxk 1   k 1 (B  k 1 A)xk  (1   k 1 )Bxk 1   k 1k 1 f
(8.4)
итерационные схемы. Здесь  k 1 ,k 1 – итерационные параметры.
Оператор B действует в том же евклидовом пространстве  n ,
что и основной оператор A :
B : u  v; u , v  n .
В дальнейшем для построения итерационных алгоритмов и
анализа их сходимости нам потребуются следующие векторы:
невязки k -й итерации
rk  A xk  f ;
(8.5)
wk  B 1rk ;
(8.6)
z k  xk  u ,
(8.7)
поправки
погрешности
где u – проекция решения краевой задачи u на сетку (в узлы) конечно-элементной модели.
Сходимость рассматривается в энергетическом подпространстве H G   n – подпространстве, порожденном в  n некоторым самосопряженным оператором G  G   0 , конструкция которого зависит от конкретного итерационного метода. В
подпространстве H G обычным образом вводятся скалярное про57
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
изведение ( x , y )G  (G x , y ) и порождаемая им норма

x
G

 G x, x  . В качестве условия окончания процесса итераций ис-
пользуется неравенство zk
G
  z0
G
, где  – заданная точность,
6
например 10 .
В схемах (8.3) и (8.4) в качестве оператора B рассмотрим линейные операторы следующих трех типов:
единичный
BE;
(8.8)
с диагональной положительно определенной матрицей
BD0;
(8.9)
B   E    0,5D  U    E    0,5 D  L   ,
(8.10)
факторизованный
где  – итерационный параметр.
Операторы D , U и L берут из представления оператора A в
виде следующей суммы:
A  L  D U ,
(8.11)
где L – оператор, соответствующий строго нижней треугольной
матрице операторы A ; D – оператор, соответствующий главной
диагонали матрицы оператора A ; U – оператор, соответствующий строго верхней треугольной матрице оператора A .
Самосопряженность и положительная определенность операторов (8.8) и (8.9) очевидны, а самосопряженность и положительная определенность факторизованного оператора (8.10) устанавливают с помощью следующей леммы.
Лемма 8.1. Если A  A  0 , то заданный выражением (8.10)
оператор B  B   0 .
Доказательство. Используя представление (8.10) оператора
B, обозначим: B1   E    0,5 D  U   и B2   E    0,5 D  L   .
Из построения следует, что B1  B2 и B2  B1 . Тогда
58
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»

B   B1B2   B2 B1  B1B2  B .
Теперь покажем, что B1  0 и B2  0 . Представим оператор A
следующим образом:
A   L  0,5 D    0,5 D  U   R1  R2 ,
где R1  L  0,5D ; R2  0,5D  U . Очевидно, что R1  R2 .
Возьмем x   n и рассмотрим скалярное произведение
 Ax, x     R1  R2  x, x    R1 x, x    R2 x, x  


  R1 x , x   x , R2 x  2  R1 x , x   2  R2 x , x   0,
откуда R1  0 и R2  0 , значит, B1  0 и B2  0 , а следовательно, и
оператор B  0 . Лемма доказана.
Известно, что операторы типа (8.8)–(8.10) ограничены [19, 40].
Определение параметров  k 1 , k 1 и  зависит от конкретного
итерационного метода. Ниже будут рассмотрены двухслойные
итерационные методы:
– метод минимальных поправок,
– метод циклических чебышевских итераций,
– метод верхней релаксации
и трехслойные итерационные методы:
– полуитерационный метод Чебышева,
– метод сопряженных поправок,
– метод сопряженных градиентов,
а также алгоритмы их реализации с разной конструкцией оператора B .
8.2. Двухслойные итерационные методы
Метод циклических чебышевских итераций. В соответствии с
общей теорией циклического чебышевского метода [19, 37] итерационные параметры k для схемы (8.3) определяют с помощью
следующей формулы:
59
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
k 
0
, k  1, N Ch ,
1  0  k
(8.12)



при  k  M N  cosi ; i  i(N ) ; i(N )  2i 1; i 1, NCh  , где M N –
2


множество корней полинома Чебышева первого рода степени N Ch ;

1 
2
0 
;   1 ; 0 
; 1 и  2 – постоянные энергетиче1 
1   2
2
ской эквивалентности операторов A и B, т. е. 1B  A   2 B ;
1  0 (в данном случае оператор G  AB 1 A ).
Для того чтобы воспользоваться формулой (8.12), необходимо
задать число итераций N Ch , определяющее длину чебышевского
цикла, и значения постоянных 1 и  2 . Число итераций в чебышевском цикле NCh можно определить из оценки нормы погрешности z N в пространстве H G , которая обычно используется в качестве условия окончания процесса итераций:
zN
где qN 
21N
1  12 N
G
 q N z0
G
  z0
, соответственно 1 
G
1 
1 
,
(8.13)
;  – требуемая точ-
ность.
Из выражения для qN имеем приближенную оценку числа
итераций:
N Ch 
ln(2 )
.
ln(1 1 )
(8.14)
После расчета числа итераций NCh в чебышевском цикле необходимо построить последовательность параметров  k , определяющую упорядоченность множества M N . Известно, что упорядоченность множества M N оказывает значительное влияние на
накопление вычислительной погрешности и, следовательно, на
60
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сходимость чебышевских итераций [19]. Упорядоченность множества M N зависит от упорядоченности множества нечетных чисел {i(N ) } . Можно использовать алгоритм построения оптимальных упорядочений множества {i( N ) } , изложенный в работе [19].
Ниже приведены примеры множеств {i( N ) } , построенных в зависимости от числа итераций N Ch :
   1, 13, 5, 9, 3, 11, 7 ;
(7)
i
   1, 15, 7, 9, 3, 13, 5, 11 ;
(8)
i
   1, 17, 7, 11, 3, 15, 5, 13, 9 ;
(9)
i
   1, 21, 9, 13, 3, 19, 7, 15, 5, 17, 11 ;
(11)
i
   1, 25, 11, 15, 5, 21, 7, 19, 3, 23, 9, 17, 13 ;
(13)
i
   1, 29, 13, 17, 5, 25, 9, 21, 3, 27, 11, 19, 7, 23, 15 ;
(15)
i
   1, 31, 15, 17, 7, 25, 9, 23, 3, 29, 13, 19, 5, 27, 11, 21 ;
(16)
i
   1, 33,15, 19, 7, 27, 9, 25, 3, 31, 13, 21, 5, 29, 11, 23, 17 ;
(17)
i
   1, 37,17, 21, 7, 31,11, 27, 3, 35,15, 23, 5, 33,13, 25, 9, 29,19 .
(19)
i
Для начальной оценки постоянных 1 и  2 энергетической эквивалентности можно использовать асимптотическое свойство
двухслойных градиентных методов [19, 37]. Выберем метод минимальных поправок, который определяется оператором G следующего вида [19, 40, 41]:
G  AB 1 A .
(8.15)
61
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для схемы (8.3) итерационные параметры k 1 в методе минимальных поправок вычисляют по формуле
k 1 
 Awk , wk 
B
1
Awk , wk

, k  0, 1, 
(8.16)
Скорость сходимости этого метода оценивается значением отношения норм векторов погрешностей двух последовательных
итераций:
k 1 
zk 1
zk
G
G

 rk 1 , wk 1 
.
 rk , wk 
(8.17)
Последовательность k  является монотонно возрастающей и
для достаточно больших k выходит на асимптоту: k    const ,
k   (в этом и состоит асимптотическое свойство). Приближенные оценки постоянных 1 и  2 спектральной энергетической эквивалентности операторов A и B могут быть найдены из квадратичного уравнения [19]
1  k  1  k 1   k k 1 .
(8.18)
Как показывает практика вычислений, для построения уравнения (8.18) достаточно выполнить четыре – семь итераций по методу минимальных поправок. Однако сам метод минимальных поправок для решения СЛАУ (8.1) по схеме (8.3) применять
нерационально из-за значительных вычислительных затрат и достаточно быстрого падения скорости сходимости итераций.
После того как из уравнения (8.18) определены значения констант 1 и  2 , по формулам (8.12) найдем параметры чебышевского итерационного цикла. Поскольку константы 1 и  2 вычисляют
приближенно, после проведения чебышевского цикла их необходимо уточнить. Для коррекции 1 и  2 можно использовать следующие соотношения [37]:
62
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 s   s Y 2  1 
1 
1s 1   1s   2s  2 1
;
2 
2
Y 
1 
 2s 1   1s   2s 
2 

2
где Y      1

2 sN
1  2s N
, s 

1
N
1  s
1  s
 2s
 1s
2
, причем  
, s 
1(s )
 (2s )
(8.19)
Y  1 
,
Y 
qNe , s
qN , s
2
,
qNe , s

zs
z0
G
, qN , s 
G
; s – номер чебышевского цикла;
N = N Ch – число итераций в чебышевском цикле; z0
G
, zs
G
–
нормы вектора погрешности после проведения чебышевского цикла с номером s .
Теоретически за NCh чебышевских итераций G -норма вектора погрешности z
уменьшается в
G
q 
e
N
должна уменьшиться в qN1 раз, реально z
1
G
раз, причем, как правило, qNe  qN . Это
происходит по двум причинам: первая – приближенность оценок
констант 1 и  2 , вторая – ошибки округления.
Практика численных исследований показывает: если после
N Ch чебышевских итераций справедливо неравенство
qN  qNe  p  1 ,
(8.20)
то перед проведением следующего чебышевского итерационного
цикла необходимо с помощью соотношений (8.19) выполнить коррекцию констант 1 и  2 . Параметр p определяется экспериментально. Проведенными численными исследованиями на разных
сетках установлено, что параметр p желательно выбирать из диапазона 0,5…0,7, причем бóльшее значение рекомендуется для
N Ch  15 [39].
63
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Если после NCh чебышевских итераций возникает ситуация,
когда qNe  1 , то необходимо выполнить одну или несколько итераций по методу минимальных поправок, чтобы подавить накопленные ошибки округления. Константы 1 и  2 в данном случае
не корректируют.
Иногда после NCh чебышевских итераций выполняется неравенство qNe  qN , которое указывает на тот факт, что практическая
скорость сходимости, определяемая qNe , не меньше теоретической,
характеризуемой значением qN . В этом случае, если решение еще
не получено, вновь повторяют чебышевский итерационный цикл с
прежними значениями констант 1 и  2 .
Результаты численных исследований показывают, что следует
проводить коррекцию только константы 1 , поскольку значение
константы  2 определяется достаточно точно из уравнения (8.18)
[37, 39].
Если в итерационной схеме (8.3) используется факторизованный оператор B, заданный соотношением (8.10), то перед проведением чебышевских итераций необходимо определить дополнительный итерационный параметр  . Между параметром  и
максимальной константой  2 существует соотношение [19]
2 
1
.
2
(8.21)
Определение констант энергетической эквивалентности 1 и
 2 операторов A и B в случае факторизованного оператора B –
более сложная задача по сравнению с ситуацией, когда оператор
B имеет единичную или диагональную матрицу. В данном случае названные константы выражаются через дополнительные
константы  и  , которые находят из неравенств, определяющих априорную информацию об операторах итерационной схемы
(8.3) [19]:
E  A ,  0,5D  U  0,5 D  L  
64

A ,   0.
4
(8.22)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Константы 1 и  2 можно выразить через  и  следующим
образом:
1 
где  


2 1 

.

При этом отношение  

, 2 

4 
,
(8.23)
1
максимально, если
2

2

.
(8.24)
Для определения констант  и  можно использовать следующий метод [19]. Положим в (8.10)   0, тогда B  E. Выполним по методу минимальных поправок N 0 итераций (обычно четыре – семь) и из уравнения
1   N 1 g1 1   N g1    N 1 N
определим константу g1  0, такую, что справедливо неравенство
2
g1E  A . Теперь положим в (8.10)    1
и снова по методу
g1
минимальных поправок выполним N 0 итераций. Затем из уравнения (8.22) найдем g 2  0 , при этом будет справедливо неравенство
g2 B  A .
Как показано в работе [19],
  g1 и   4
g1  g 2  g1 g 2 1
g1 g 2 12
.
(8.25)
Значение  вычисляется по  и  с помощью выражения
(8.24). Если справедливо неравенство
2 g 2  g1 (1  g 2  1) ,
(8.26)
65
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
то значение  будет принадлежать интервалу  0, 1  .
Когда неравенство (8.26) не выполняется, необходимо несколько увеличить значение параметра 1 и провести все расчеты
заново (рекомендуется брать 1 2 g 2 ).
Таким образом, решение СЛАУ (8.1) осуществляется комбинированным методом, сочетающим итерации, выполняемые методом
минимальных поправок, и чебышевские итерационные циклы.
Алгоритм рассмотренного комбинированного метода состоит
из следующих основных вычислительных шагов:
1) если матрица оператора B является диагональной (8.8) или
(8.9), то проводим N 0 итераций ( N 0 = 4…7) по методу минимальных поправок и из уравнения (8.18) определяем константы 1 и
 2 ; если оператор B имеет факторизованный вид (8.10), то проводим две серии из N 0 итераций по методу минимальных поправок,
из (8.25) определяем величины  и  , а затем из (8.23) и (8.24)
находим значения констант 1 ,  2 и параметра  ;
2) с помощью соотношений (8.12) вычисляем итерационные
параметры k чебышевского цикла, а затем проводим цикл из
N Ch чебышевских итераций и проводим оценку достигнутой точности;
3) если требуемая точность не достигнута, по результатам чебышевского цикла оцениваем значение qNe ;
4) выполняем сравнение qNe и qN :
при p  qNe (или 1  qNe ) проводим N MP итераций методом минимальных поправок, оцениваем достигнутую точность и, если
требуемая точность не достигнута, выполняем п. 2 алгоритма;
если qNe  qN , то вновь выполняем п. 2 алгоритма;
если qN  qNe  p  1 , то в соответствии с (8.19) проводим коррекцию константы 1 (одновременно корректируем константу  2
и параметр  ), а затем выполняем п. 2 алгоритма;
Численные исследования показывают, что комбинированный
метод, сочетающий чебышевские итерационные циклы и итерации
по методу минимальных поправок, обладает существенно более вы66
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сокой сходимостью по сравнению с решением только по методу
минимальных поправок. На скорость сходимости комбинированного метода большое влияние оказывает длина чебышевского цикла
N Ch . Число NCh , как отмечалось выше, может быть оценено по
формуле (8.18), однако на практике число чебышевских итераций в
цикле выбирают не очень большим ( N Ch  7…17). Это связано с
тем, что накопление ошибок округления при больших N Ch могут
вызвать затруднения с обеспечением численной устойчивости.
Большое влияние на сходимость оказывает конструкция оператора B. Использование единичного оператора B в схемах (8.3) и
(8.4) для данных методов (чебышевские итерационные циклы и
метод минимальных поправок) делает расчеты практически невозможными – крайне низкая сходимость. При одной и той же длине
чебышевского цикла NCh факторизованный оператор B обеспечивает существенно более высокую скорость сходимости по общему числу итераций. Преимущества факторизованного оператора
по времени решения задачи заметны только с ростом числа узлов
сеток конечно-элементных моделей.
Метод верхней релаксации. Метод реализуется двухслойной
итерационной схемой (8.3), если в ней положить
B  D  L ,  k   .
(8.27)
Операторы D и L берут из (3.11). Оптимальное значение итерационного параметра  определяется спектральным радиусом оператора перехода S  E    D  L 
1
A , однако на практике 
подбирают экспериментально из диапазона  0, 2  .
Метод верхней релаксации на всех сетках и по общему числу
итераций N и по времени счета  превосходит комбинированный
метод с диагональным оператором B  D. По числу итераций N
он уступает комбинированному методу с факторизованным оператором B, но по времени счета вполне сопоставим на малых и
средних сетках. На больших сетках преимущества комбинированного метода достигаются только при чебышевских циклах с числом итераций NCh = 13, 15.
Таким образом, из-за простоты реализации метод верхней релаксации в классе двухслойных схем вполне применим для ре67
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
шения СЛАУ с разреженной матрицей с числом неизвестных до
5000 – 7000.
8.3. Трехслойные итерационные методы
Если оператор B  E , то неявная стандартная трехслойная
схема имеет вид, заданный выражением (8.4), причем первую итерацию выполняют по двухслойной схеме (8.3). Конкретные формулы для итерационных параметров k и  k определяет тот или
иной итерационный метод. Ниже рассмотрены некоторые трехслойные методы. В качестве оператора B использованы операторы, заданные соотношениями (8.9) и (8.10).
Полуитерационный метод Чебышева. В этом методе итерационные параметры определяют следующим образом [19, 37]:
k 
2
4
,  k 1 
, k = 1, 2, ...; 1 = 2,
1   2
4  02  k
(8.28)

1 
, причем   1 .
1 
2
Константы 1 и  2 берут из соотношения энергетической эквивалентности операторов A и B : 1B  A   2 B , где 1  0 (в
где 0 
данном случае оператор G  AB 1 A ). Для численного определения
1 и  2 используют двухслойную схему (8.3) и метод минимальных поправок, обладающий асимптотическим свойством, что позволяет применить соотношение (8.18).
Данный метод не является циклическим, и число итераций
NCh может быть любым, при этом параметр k не меняется.
Однако на практике проводят некоторое заданное число итераций NCh , а затем выполняют анализ сходимости. При необходимости делают коррекцию констант 1 и  2 . Алгоритм решения СЛАУ (8.1) в данном случае полностью совпадает с тем,
который был рассмотрен в разд. 8.2 применительно к чебышевскому методу (8.12). Отличие касается вычисления итерационных параметров – п. 2 алгоритма. Вычисление необходимо проводить по формулам (8.28).
68
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для фиксированного числа итераций NCh оба метода – полуитерационный Чебышева и чебышевский двухслойный – имеют
одинаковую оценку уменьшения нормы вектора погрешности zk ,
заданного соотношением (8.13). Однако эта оценка для двухслойного метода верна, если выполнены все NCh итераций, а для трехслойного метода эта оценка справедлива и для всех промежуточных итераций. На промежуточных итерациях у трехслойного
метода наблюдается монотонное уменьшение нормы вектора погрешности в отличие от двухслойного, у которого могут происходить некоторые колебания нормы вектора погрешности. Таким
образом, вычислительная устойчивость полуитерационного метода
Чебышева существенно выше, чем у чебышевского двухслойного
циклического метода.
При длинных циклах ( N Ch достаточно большое, например 15
и более итераций) общее время решения может оказаться меньше,
чем при коротких циклах, хотя общее число итераций N для достижения заданной точности, больше. Это зависит от необходимости подавления ошибок округления с помощью итераций по методу минимальных поправок. Итерации подавления ошибок
округления методом минимальных поправок используют, когда
справедливы неравенства (8.20) (см. п. 4 алгоритма двухслойного
метода).
Очевидно, что за большое число итераций N Ch (с учетом монотонности сходимости полуитерационного метода) значение qNe
становится меньше значения порога p. За короткие циклы этого
не происходит, поэтому включают итерации подавления, которые
в силу конструкции метода минимальных поправок выполняются
значительно дольше чебышевских. Это явление наблюдается как
для оператора B с диагональной матрицей, так и для факторизованного оператора B. Использование единичного оператора B
(8.8) в полуитерационном методе Чебышева является крайне неэффективным.
Как показывают численные исследования, по числу итераций
оператор B с диагональной матрицей существенно уступает факторизованному. По времени счета преимущество также у факторизованного оператора. Полуитерационный метод Чебышева с факторизованным оператором по числу итераций имеет заметное
69
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
преимущество перед методом верхней релаксации, однако по времени счета этот метод более эффективен только на достаточно
больших сетках (более 7 000–10 000 узлов) при проведении циклов
длиной N Ch  8…17 итераций с коррекцией констант 1 и  2 и
подавлением ошибок округления.
Метод сопряженных поправок. Метод – частный случай более
общего – метода сопряженных направлений, который, в свою очередь, является одним из методов вариационного типа. В данном
методе итерационные приближения проводятся по общей трехслойной схеме (8.4). К операторам схемы предъявляют следующие
требования [19, 40, 41]:
A  A  0; B  B  0;
(8.29)
1B  A   2 B ; 1  0.
Энергетическое пространство H G вводится с помощью оператора G  AB 1 A . Итерационные параметры k 1 и  k 1 определяют по формулам [19]
k 1 
 Awk , wk 
B
1
Awk , wk

;
1
(8.30)
 
(Awk , wk ) 1 
 k 1  1  k 1
 .
k (Awk 1 , wk 1 )  k 

Для метода справедлива оценка (8.13), указывающая на то, что
метод сопряженных поправок, вообще говоря, сходится не хуже
двухслойного чебышевского, но при этом имеет место монотонность убывания нормы вектора погрешности zk G . При рассмотрении двухслойных итерационных схем (8.3) методу сопряженных
поправок соответствует метод минимальных поправок. Применение
трехслойных итерационных схем, особенно неявных, оправдывается
тем, что они обеспечивают более высокую скорость сходимости.
При одной и той же конструкции оператора B трехслойные
схемы для метода сопряженных поправок обладают значительно
70
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
более высокой сходимостью и по числу итераций N и по времени
счета  на сетках с малым и средним числом узлов (менее 5 000 –
7 000). Однако независимо от конструкции оператора B метод
сопряженных поправок для сеток с достаточно большим числом
узлов (более 40 000) оказывается совершенно неэффективным изза крайне низкой сходимости.
Метод сопряженных градиентов. В настоящее время метод
сопряженных градиентов – один из самых перспективных итерационных методов для решения СЛАУ больших порядков, в том
числе и с плохо обусловленными матрицами [19, 37, 42]. Кроме
того, он может быть использован для ускорения других итерационных методов [37]. Для реализации метода сопряженных градиентов требуется выполнение условий (8.29), а энергетическое пространство H G вводится оператором G  A . Тогда итерационные
приближения вычисляют по схеме (8.4), а итерационные параметры k 1 и  k 1 находят по следующим формулам [37, 38]:
k 1 
 rk , wk 
 Awk , wk 
, k  0, 1, .... ;
 
 rk , wk  1
 k 1  1  k 1

k  rk 1 , wk 1   k

1

 , k  1, 2, ...; 1  1.

(8.31)
Для метода сопряженных градиентов сохраняется оценка
(8.13). Алгоритм его реализации совпадает с алгоритмом для метода сопряженных поправок, рассмотренным выше. В случае
применения двухслойной схемы (8.3) методу сопряженных градиентов соответствует метод скорейшего спуска – итерационные
параметры k у этих методов вычисляются по одной формуле.
Кроме того, если B  B , то метод скорейшего спуска обладает
асимптотическим свойством и его можно использовать при факторизованном операторе B для определения параметра  .
Как и во всех предыдущих случаях, усложнение конструкции
оператора B ведет к сокращению общего числа итераций. Особенно заметное уменьшение числа итераций дает применение
факторизованного оператора B на средних и больших сетках.
71
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Сокращение числа итераций при использовании метода сопряженных градиентов ведет и к уменьшению времени счета, за исключением сеток с малым числом узлов. Это связано с тем, что
процедуры метода сопряженных градиентов при факторизованном операторе B выполняются дольше по сравнению с аналогичными, использующими операторы B с диагональными матрицами. Таким образом, заметный выигрыш во времени
происходит, только если удается добиться значительного сокращения общего числа итераций, как, например, это происходит на
средних и больших сетках.
По сравнению со всеми другими метод сопряженных градиентов на всех сетках и для любых конструкций оператора B
имеет явное преимущество по числу итераций и по времени их
реализации.
8.4. Локально оптимальные трехслойные методы
Формулы для вычисления итерационных параметров k 1 и
 k 1 рассмотренных выше методов сопряженных поправок и сопряженных градиентов были построены для идеального вычислительного процесса, поэтому они не учитывают в достаточной мере
влияние ошибок округления [19]. Формулы (8.30) и (8.31) получены
исходя из минимизации нормы вектора погрешности zk 1 G в пространстве H G при переходе от x0 к xk 1 . Существует и другой
подход к построению итерационных параметров k 1 и  k 1 . Пусть
известны приближенные решения xk 1 и xk , тогда формулы для
итерационных параметров можно построить из условия достижения
минимума нормы вектора погрешности zk 1 G в пространстве H G
за один шаг по трехслойной схеме, т. е. из условия локальной оптимизации при переходе от k-й к (k + 1)-й итерации. Как показали численные исследования, такой подход обеспечивает существенно
бóльшую численную устойчивость итерационных методов вариационного типа и соответственно быструю сходимость.
Итерационные параметры k 1 и  k 1 для локально оптимальных трехслойных методов вариационного типа вычисляют по
формулам [19]
72
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
 k 1 
 ak  bk  bk  d k ek
,
 ck  d k  ek   ak  bk 2
k  1, 2,  ; 1  1;
(8.32)
a 1   k 1 bk
k 1  k 
, k  0 , 1, 2, 
ek
 k 1 ek
Здесь приняты следующие обозначения для скалярных произведений:
ak   Gwk , zk  ;
bk   Gwk , zk 1  ;
ck   Gzk , xk  xk 1  ;
(8.33)
d k   Gzk 1 , xk  xk 1  ;
ek   Gwk , wk  .
Выбор конкретного метода из всего семейства локально оптимальных зависит от структуры оператора G .
Метод сопряженных поправок. Для реализации этого метода
в рамках трехслойной локально оптимальной схемы операторы A
и B должны удовлетворять условиям (8.29). Энергетическое пространство H G вводится с помощью оператора G  AB 1 A . Тогда
скалярные произведения (8.33) примут следующий вид:
ak   Awk , wk  ;
bk   Awk , wk 1  ;
ck   wk , rk  rk 1  ;
(8.34)
d k   wk 1 , rk  rk 1  ;


ek  B 1 Awk , Awk .
Итерационные параметры k 1 и  k 1 вычисляют по формулам (8.32).
73
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Численные исследования показывают, что локально оптимальная схема обеспечивает сходимость за существенно меньшее число итераций и соответственно меньшее время счета. Особенно заметно преимущество перед двухслойной схемой. Кроме того,
локально оптимальная схема позволяет получить решение даже на
достаточно больших сетках за малое число итераций и вполне
приемлемое время. Двухслойная и обычная трехслойная схемы для
методов соответственно минимальных и сопряженных поправок
этого сделать не позволяют.
При реализации локально оптимальной схемы факторизованный оператор B обеспечивает сходимость за меньшее число итераций по сравнению с оператором B  D , имеющим диагональную матрицу. Однако применение факторизованного оператора B
не дает выигрыша во времени решения СЛАУ. Очевидно, это связано с большими вычислительными затратами на каждую итерацию при использовании факторизованного оператора.
Метод сопряженных градиентов. Использование локально
оптимальной схемы для метода сопряженных градиентов предполагает, что операторы A и B удовлетворяют требованиям (8.29), а
энергетическое пространство H G вводится оператором G  A .
Тогда скалярные произведения (8.33) примут вид
ak  (wk , rk );
bk  (wk , rk 1 );
ck  (rk , xk  xk 1 );
(8.35)
d k  (rk 1 , xk  xk 1 );
ek  (Awk , wk ) .
Для вычисления итерационных параметров k 1 и  k 1 используют формулы (8.32).
В этой схеме, как и ранее, усложнение конструкции оператора
B ведет к уменьшению числа итераций и времени решения СЛАУ.
Исключение составляют, как всегда, сетки с малым числом узлов,
где уменьшение числа итераций не дает выигрыша во времени при
использовании факторизованного оператора B .
74
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Сравнение метода сопряженных градиентов, реализованного
в трехслойной локально оптимальной схеме, со всеми рассмотренными выше методами в двух- и трехслойных схемах, указывает на значительное его преимущество и по общему числу итераций N , и по времени решения СЛАУ  . Существенное влияние
на число итераций N и время решения  оказывает конструкция
оператора B . На сетках, имеющих более 5000 узлов, использование факторизованной конструкции оператора B представляется
вполне оправданным. При решении задач с меньшим числом узлов может быть использован оператор B с диагональной матрицей. Сравнительный анализ вычислительной эффективности прикладных итерационных методов решения сеточных уравнений
сделан в работе [42].
8.5. Построение и использование разреженных матриц
Матрицы СЛАУ, возникающие в процессе построения дискретных аналогов краевых задач с помощью метода конечных
элементов, как правило, симметричны с выраженной разреженной
структурой. При соответствующей нумерации узлов конечноэлементных моделей эти матрицы могут иметь ленточную структуру. В связи с этим при использовании конечно-элементной технологии возникает проблема разработки эффективных алгоритмов
формирования, хранения и применения разреженных матриц.
Этим вопросам посвящена обширная литература [22, 26–37].
Память, используемая для хранения разреженных матриц, состоит из двух частей: основной, в которой содержатся числовые
значения элементов матриц, и дополнительной, где хранятся указатели, индексы и другая информация, необходимая для формирования структуры матриц и обеспечивающая доступ к численным
значениям элементов матриц при формировании и решении
СЛАУ, т. е. так называемые списки связности. Способы хранения
и использования данных, размещаемых в основной и дополнительной памяти, весьма разнообразны и зависят, главным образом,
от метода решения СЛАУ. Для структурного построения информации, хранящейся в дополнительной памяти, – списков связности – широко используют теорию графов [31]. Особенно сложные
алгоритмы работы с разреженными матрицами возникают при
прямых методах решения, например типа фронтальных. В первую
75
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
очередь это связано с тем, что при разложении разреженной матрицы A  LLT она обычно претерпевает некоторое заполнение, т. е.
нижний треугольный множитель L имеет ненулевые элементы в
тех позициях, где в исходной матрице A стояли нули; не исключено, что при этом могут появиться и новые нулевые элементы.
Подобных проблем при использовании итерационных методов не
возникает. Кроме того, прямые методы требуют такого упорядочения нумерации узлов конечно-элементной сетки, которое обеспечивает минимальную ширину ленты.
Ниже приведен алгоритм построения и использования разреженных матриц, отличающийся простотой и высокой вычислительной эффективностью [22, 39].
Рассмотрим в качестве примера конечно-элементную модель, представленную на рис. 8.1. С каждым узлом сетки связано некоторое количество узлов, которые могут быть разбиты на два
Рис. 8.1. Конечно-элементная модель
множества. Первое множество составляют узлы, номера которых
меньше i – номера рассматриваемого узла, а второе множество –
узлы, номера которых больше i . На рис. 8.2 показана структура
матрицы СЛАУ, соответствующая узловым связям конечноэлементной модели (см. рис. 8.1). Матрица A хранится в виде
двух векторов (рис. 8.3): первый вектор d – члены главной диагонали, второй вектор u – ненулевые элементы, расположенные
над главной диагональю. Для формирования и использования
элементов второго вектора создают дополнительные четыре вектора, содержащие информацию о связях узлов конечно-элементной модели.
76
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1
2
3
1
a11
2
a 21
a 31
a12
a 22
a 32
a13
a 23
a 33
3
a44
a54
a 64
a 74
a84
a 94
4
5
6
a51
a 61
4
a52
a 62
7
8
9
5
6
a15
a 25
a16
a 26
a 45
a55
a 65
a 46
a56
a 66
a 95
a 96
Рис. 8.2. Структура матрицы
1
2
3
4
5
6
7
8
9
d
a11
a 22
a 33
a 44
a 55
a 66
a 77
a 88
a 99
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
7
8
9
a 47
a 48
a 77
a87
a 97
a 78
a88
a 98
a 49
a59
a 69
a 79
a89
a 99
A
u
a15
a16
a12
a13
a 25
a 26
a 23
a 47
a 48
a 49
a 46
a 45
a 59
a 56
a 69
a 78
a 79
a89
Рис. 8.3. Структура векторного хранения матрицы A :
d – вектор элементов главной диагонали; u – вектор наддиагональных элементов
77
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рассмотрим i-е уравнение СЛАУ (8.1). В общем виде оно имеет
вид
i 1
 aij x j  aii xi 
j 1
N
 aij x j  fi ,
(8.36)
j i 1
где aij – элементы матрицы A ; fi – элемент вектора правой части
f ; x j – элементы вектора неизвестных x ; N – число неизвестных.
В уравнении (8.36) содержатся суммы произведений двух типов. Первая сумма состоит из произведений элементов aij матрицы A и элементов x j вектора неизвестных x , у которых j  i , а
для произведений второй суммы имеем j  i . Для построения произведений первой суммы используют два вектора: v NND и v NUD ,
пример структуры которых показан на рис. 8.4. В векторе v NND
хранятся номера тех узлов, связанных с рассматриваемым узлом, у
которых номера меньше номера рассматриваемого узла. Например, рассмотрим узел с номером 5 (см. рис. 8.1). Этот узел входит
в элементы с номерами 2 и 3. В силу этого он имеет связь с узлами
4, 2 и 1, номера которых меньше 5. Эти узлы записывают в вектор
v NUD в последовательности, определенной номерами элементов и
правилом обхода их узлов. На рис. 8.1 стрелками показано начало
и направление обхода узлов элементов. Этим объясняется порядок
появления узлов 4, 2 и 1 в векторе v NUD . В векторе v NND хранятся
номера индексов, координирующих в векторе v NUD расположение
первого узла всего массива номеров узлов, связанных с рассматриваемым. Например, как видно из рисунка, массив узлов связанных
с узлом 5, начинается с индекса 6. Последний индекс массива определяется начальным индексом узла, следующего за рассматриваемым, за вычетом единицы, например, для узла 5 последний индекс массива в векторе v NUD равен 8 = 9 – 1, где 9 – начало
массива узлов для узла 6.
Аналогично строят еще два вектора v NNP и vNUP , содержащих
информацию о номерах узлов, связанных с рассматриваемым, но
номера которых больше номера рассматриваемого. На рис. 8.4
78
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
показан пример конструкции указанных векторов применительно
к конечно-элементной модели, представленной на рис. 8.1. В векторе v NNP указывается начальный индекс, фиксирующий начало
списка узлов, размещенных в векторе vNUP и связанных с рассматриваемым.
v NND
1
2
3
4
5
6
7
8
9
18
1
2
0
4
7
11
12
14
v NUD
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
1
1
2
4
2
1
4
5
2
1
4
4
7
7
8
4
6
5
v NNP
1
2
3
4
5
6
7
8
1
5
0
8
13
15
1
18
v NUP
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
5
6
2
3
5
6
3
7
8
9
6
5
9
6
9
8
9
9
Рис. 8.4. Структура векторов-указателей
В тех редких случаях, когда рассматриваемый узел конечноэлементной модели не имеет связи с узлами, номера которых
меньше или больше номера рассматриваемого узла, в соответствующие ячейки векторов v NND или v NNP заносят нули. Например,
узел 4 (см. рис. 8.1) не имеет связи с узлами, номера которых меньше 4, это учитывается записью нуля в ячейку 4 вектора v NND (см.
рис. 8.4). Узел 3 не имеет связи с узлами, номера которых больше
3, этому обстоятельству соответствует запись нуля в ячейку 3 вектора v NNP . Исключение составляет запись в ячейку 1 вектора
v NND , куда заносится номер последней ячейки вектора v NUD . Это
79
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сделано для того, чтобы зафиксировать параметры оператора цикла при просмотре списка узловых связей последнего узла сетки
конечно-элементной модели. Для всех других узлов эти параметры
определяются тривиально. Для каждой конечно-элементной модели векторы строят один раз – сразу по окончании построения массива, описывающего конечные элементы сетки.
Структура вектора u (см. рис. 8.3) в значительной степени повторяет структуру вектора vNUP . Ненулевые наддиагональные элементы aij матрицы A хранятся в векторе u построчно, а порядок
их следования в пределах строки определяется порядком следования индексов в векторе vNUP .
Теперь рассмотрим построение сумм произведений, входящих
в уравнение (8.36). Начнем со второй суммы:
S Pi 
N
 aij x j  aii 1xi 1  aii 2 xi 2  ...  aiN xN .
(8.37)
j i 1
При фиксированном индексе i определяют номера первой и
последней ячеек вектора v NNP , где размещаются номера узлов,
связанных с i-м, для которых j  i . Например, для узла 5 (см.
рис. 8.1) имеем: номер первой ячейки – 13, а номер последней
ячейки 14 = 15 – 1. Номера первой и последней ячеек определяют
в векторе u диапазон размещения ненулевых элементов матрицы
A , участвующих в произведениях суммы (8.37). Кроме того, эти
же номера ячеек определяют в векторе vNUP диапазон размещения
номеров соответствующих элементов xi вектора x . Таким образом, например, для узла 5 сумма (8.37) состоит из двух слагаемых
и имеет вид
S P5  a59 x9  a56 x6 .
(8.38)
Построение первой суммы из уравнения (8.36)
i 1
S Di   aij x j  ai1 x1  ai 2 x2  ...  aii 1 xi 1
j 1
80
(8.39)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
алгоритмически сложнее. При фиксированном номере i определяют номера первой n и последней m  n  m  ячеек вектора
v NND , где хранятся номера узлов j , связанных с узлом i , для которых j  i . Этим самым задаются компоненты xi вектора x . Теперь необходимо определить элементы aij  a ji
i  j 
матрицы
A , хранящиеся в векторе u . Эта процедура выполняется следующим образом. Последовательно из ячеек с n по m в векторе v NUD
 j  i  . По этим номерам j из вектора v NNP
определяют номера первой k и последней l  l  k  ячеек вектора
v NUP . Для каждого узла j в диапазоне ячеек  k  l  вектора v NUP
берут номера узлов j
находят номер i , поскольку, во-первых, узел i связан с каждым
узлом j , а, во-вторых, i  j . Просматривая ячейки вектора v NUP и
сравнивая номера узлов, находящихся в этих ячейках, с номером
i , можно найти номер ячейки, в которой для данного узла j и
диапазона номеров его ячеек  k  l  находится номер узла i . По
номеру ячейки узла i в векторе v NNP можно найти в векторе u
элемент aij  a ji
i 
j , j  i  , т. е. точно так же, как это делалось
при построении второй суммы (8.36).
Пример. Построим сумму (8.39) для узла 5 конечно-элементной модели, показанной на рис. 8.1. С этим узлом связаны узлы
4, 2 и 1. Номера их хранятся в векторе v NUD в ячейках с 4-й по
6-ю, причем они меньше 5. Таким образом, компоненты xi вектора x известны: x4 , x2 и x1 . Определим элементы матрицы
A . Узлу 4 в векторе v NUP соответствуют ячейки с 8-й по 12-ю.
В этих ячейках хранятся номера узлов, связанных с номером 4,
причем все эти номера больше 4. Один из этих номеров 5. Этот
номер (узла) находится в ячейке вектора v NUP с номером 12. В
ячейке с этим же номером в векторе u хранится элемент
a45  a54  j  4  i  5  матрицы A . Аналогично отыскиваются
элементы a25  j  2  i  5  и a15
получаем выражение
 j  1  i  5 .
Таким образом,
81
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
S D5  a45 x4  a25 x2  a15 x1 .
(8.40)
Окончательно, с учетом (3.38), вся сумма S 5 приобретет следующий вид:
S 5  S D5  a55 x5  S P5  a45 x4  a25 x2 
 a15 x1  a55 x5  a59 x9  a56 x6 .
(8.41)
Из представленного алгоритма следует, что вычислительные
затраты на построение сумм произведений S Di и S Pi разные. Вторая сумма S Pi строится значительно быстрее. Это связано с тем,
что при построении первой суммы S Di необходимо проводить дополнительный сравнительный анализ номеров узлов, размещенных
в векторе v NUP . При рассмотрении процедур алгоритма предполагалось, что все необходимые данные располагаются в оперативной
памяти компьютера. Длина одномерных массивов, в которых размещаются векторы u , v NND и v NNP , зависит от структуры рассматриваемой конечно-элементной модели.
Изложенный выше алгоритм формирования сумм в уравнении
(8.36) может быть применен и для построения произведений типа
Aw и B 1 Aw , используемых в ранее рассмотренных методах решения СЛАУ.
82
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ЛИТЕРАТУРА
1. Власова Е.А. Приближенные методы математической физики:
Учеб. для вузов / Власова Е.А., Зарубин В.С., Кувыркин Г.Н.; под ред.
В.С. Зарубина, А.П. Крищенко. М.: Изд-во МГТУ им. Н.Э. Баумана, 2001.
(Сер. Математика в техническом университете; Вып. XIII).
2. Зенкевич О. Метод конечных элементов в технике / Зенкевич О. М.:
Мир, 1975.
3. Сегерлинд Л. Применение метода конечных элементов / Сегерлинд Л. М.: Мир, 1979.
4. Шабров Н.Н. Метод конечных элементов в расчетах тепловых двигателей / Шабров Н.Н. Л.: Машиностроение, 1983.
5. Метод конечных элементов в механике твердых тел / А. Сахаров,
В. Кислоокий, В. Киричевский, И. Альтенбах, У. Габберт, Ю. Данкерт,
Х. Кеплер, З. Кочык; под ред. А. Сахарова, И. Альтенбаха. Киев: Вища
школа, 1982.
6. Зенкевич О. Конечные элементы и аппроксимация / Зенкевич О.,
Морган К. М.: Мир, 1986.
7. Зарубин В.С. Инженерные методы решения задач теплопроводности / Зарубин В.С. М.: Энергоатомиздат, 1983.
8. Зарубин В.С., Станкевич И.В. Расчет теплонапряженных конструкций / Зарубин В.С., Станкевич И.В. М.: Машиностроение, 2005.
9. Зарубин В.С., Математические модели механики и электродинамики сплошной среды / Зарубин В.С., Кувыркин Г.Н. М.: Изд-во МГТУ им.
Н.Э. Баумана, 2008.
10. Мартинсон Л.К., Малов Ю.И. Дифференциальные уравнения математической физики: Учеб. для вузов / Мартинсон Л.К., Малов Ю.И. под
ред. В.С. Зарубина, А.П. Крищенко. М.: Изд-во МГТУ им. Н.Э. Баумана,
1996. 363 с. (Сер. Математика в техническом университете; Вып. XII).
11. Кошелев А.И. О сходимости метода последовательных приближений для квазилинейных эллиптических уравнений // Докл. АН СССР.
1962. Т. 142, № 5. С. 1007–1010.
12. Станкевич И.В. Сходимость метода простых итераций при решении нелинейных стационарных уравнений теплопроводности //
Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. 1995. № 3.
С. 97 – 102.
83
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
13. Свешников А.Г. Лекции по математической физике / Свешников А.Г.,
Боголюбов А.Н., Кравцов В.В. М.: Изд-во Московского университета,
1993.
14. Ладыженская О.А. Краевые задачи математической физики / Ладыженская О.А. М.: Наука, 1973.
15. Стренг Г. Теория метода конечных элементов / Стренг Г., Фикс Дж.
М.: Мир, 1977.
16. Сьярле Ф. Метод конечных элементов для эллиптических задач /
Сьярле Ф. М.: Мир, 1980.
17. Епифанов В.М., Станкевич И.В. Математическое моделирование
трехмерных температурных полей турбинных лопаток // Докл. РАН.
1993. Т. 330, № 3. – С. 318–320.
18. Епифанов В.М., Станкевич И.В. Математическое моделирование
температурных полей сопловых и рабочих лопаток газотурбинных двигателей // Вестник МГТУ им. Н.Э. Баумана. 1993. № 1. С. 46–55.
19. Самарский А.А. Методы решения сеточных уравнений / Самарский А.А., Николаев Е.С. М.: Наука, 1978.
20. Ректорис К. Вариационные методы в математической физике и
технике / Ректорис К. М.: Мир, 1985.
21. Марчук Г.И. Введение в проекционно-сеточные методы / Марчук Г.И., Агошков В.И. М.: Наука, 1981.
22. Станкевич И.В. Хранение и использование разреженных матриц в
конечно-элементной технологии // Информационные технологии. 1998.
№ 12. С. 9–12.
23. Бахвалов Н.С. Численные методы / Бахвалов Н.С., Жидков Н.П.,
Кобельков Г.М. М.: Наука, 1987.
24. Кувыркин Г.Н. Термомеханика деформируемого твердого тела
при высокоинтенсивном нагружении / Кувыркин Г.Н. М.: Изд-во МГТУ
им. Н.Э. Баумана, 1993.
25. Уманский С.Э. Оптимизация приближенных методов решения
краевых задач механики / Уманский С.Э. Киев: Наукова думка, 1983.
26. Зарубин В.С. Математическое моделирование в технике: Учеб.
для вузов / под ред. В.С. Зарубина, А.П. Крищенко. М.: Изд-во МГТУ им.
Н.Э. Баумана, 2001. 496 с. (Сер. Математика в техническом университете;
Вып. XXI, заключительный).
27. Эстербю О. Прямые методы для разреженных матриц / Эстербю О.,
Златев З. М.: Мир, 1987.
28. Тьюарсон Р. Разреженные матрицы / Тьюарсон Р. М.: Мир, 1977.
29. Брамеллер А. Слабозаполненные матрицы / Брамеллер А., Аллан Р.,
Хэмэм Я. М.: Энергия, 1979.
30. Райс Дж. Матричные вычисления и математическое обеспечение /
Райс Дж. М.: Мир, 1984.
31. Джордж А. Численное решение больших разреженных систем
уравнений / Джордж А., Лю Дж. М.: Мир, 1984.
84
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
32. Икрамов Х.Д. Численные методы для симметричных линейных
систем / Икрамов Х.Д. М.: Наука, 1988.
33. Ильин В.П. Методы неполной факторизации для решения алгебраических систем / Ильин В.П. М.: Физматлит, 1995.
34. Писсанецки С. Технология разреженных матриц / Писсанецки С.
М.: Мир, 1988.
35. Парлетт Б. Симметричная проблема собственных значений.
Численные методы / Парлетт Б. М.: Мир, 1983.
36. Уилкинсон Дж. Х., Райнш К. Справочник алгоритмов на языке
АЛГОЛ. Линейная алгебра / Уилкинсон Дж. Х., Райнш К. М.: Машиностроение, 1976.
37. Хейгеман Л., Янг Д. Прикладные итерационные методы / Хейгеман Л., Янг Д. М.: Мир, 1986.
38. Дьяконов Е.Г. Минимизация вычислительной работы. Асимптотически оптимальные алгоритмы для эллиптических задач / Дьяконов Е.Г.
М.: Наука, 1989.
39. Станкевич И.В. Сравнительный анализ вычислительной эффективности прикладных итерационных методов решения сеточных уравнений теплопроводности // Тепломассообмен: Труды второй Российской
национальной конференции. М., 1998. Т. 7. С. 213–216.
40. Воеводин В.В. Матрицы и вычисления / Воеводин В.В., Кузнецов Ю.А. М.: Наука, 1984.
41. Малышев А.Н. Введение в вычислительную линейную алгебру /
Малышев А.Н. Новосибирск: Наука, 1991.
42. Некоторые современные методы решения сеточных уравнений /
А.А. Самарский, И.Е. Капорин, А.Б. Кучеров, Е.С. Николаев // Изв. вузов.
Сер. Математика. 1983. № 7. С. 3–12.
85
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ОГЛАВЛЕНИЕ
Введение ......................................................................................................
1. Постановка нелинейной стационарной задачи теплопроводности .......
2. Вариационная формулировка стационарной задачи теплопроводности .................................................................................................................
3. Постановка нелинейной нестационарной задачи теплопроводности .....
4. Построение матричных соотношений МКЭ .........................................
4.1. Понятие конечного элемента .......................................................
4.2. Стационарная задача теплопроводности .....................................
4.3. Нестационарная задача теплопроводности .................................
5. Изопараметрические отображения и функции формы конечных элементов ..........................................................................................................
5.1. Построение изопараметрических отображений .........................
5.2. Функции формы конечных элементов ........................................
6. Особенности численного интегрирования матричных соотношений
МКЭ ..............................................................................................................
6.1. Интегрирование по объему ...........................................................
6.2. Интегрирование по поверхности .................................................
7. Особенности численного решения задачи Коши .................................
7.1. Двухслойные разностные схемы ..................................................
7.2. Трехслойные разностные схемы ..................................................
7.3. Диагонализация матрицы теплоемкости .....................................
8. Методы решения систем линейных алгебраических уравнений ........
8.1. Основные понятия теории итерационных методов ....................
8.2. Двухслойные итерационные методы ...........................................
8.3. Трехслойные итерационные методы ...........................................
8.4. Локально оптимальные трехслойные методы ............................
8.5. Построение и использование разреженных матриц ...................
Литература ...................................................................................................
86
3
5
7
10
12
12
18
24
27
27
30
36
37
43
47
47
48
53
55
55
59
68
72
75
83
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Учебное издание
Котов Александр Валерианович
Станкевич Игорь Васильевич
РЕШЕНИЕ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ
МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
Редактор В.М. Царев
Корректор Е.В. Авалова
Компьютерная верстка А.Ю. Ураловой
Подписано в печать 14.12.2010. Формат 60×84/16.
Усл. печ. л. 5,12. Тираж 600 экз. Изд. № 3. Заказ
Издательство МГТУ им. Н.Э. Баумана.
Типография МГТУ им. Н.Э. Баумана.
105005, Москва, 2-я Бауманская ул., 5.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Для заметок
Документ
Категория
Другое
Просмотров
565
Размер файла
730 Кб
Теги
решение, методов, конечный, элементов, теплопроводность, задачи
1/--страниц
Пожаловаться на содержимое документа