Использование пакета конечных элементов FreeFen++ для задач гидродинамики, электрофореза и биологии

Южный федеральный университет

М. Ю. Жуков, Е. В. Ширяева

Использование пакета
конечных элементов FreeFem++
для задач гидродинамики,
электрофореза и биологии

Ростов-на-Дону

2008

УДК 519.63; 532.5; 537.364
ББК 22.25

Ж 86

Печатается по решению
кафедры вычислительной математики и математической физики

и Совета факультета математики, механики и компьютерных наук ЮФУ
от 15.09.2007

Р е ц е н з е н т ы

Я. М. Ерусалимский, Н. В. Петровская

Жуков М. Ю., Ширяева Е. В.

Использование пакета конечных элементов FreeFem++ для задач гидродинамики, электрофореза и биологии. — Ростов н/Д :
Изд-во ЮФУ, 2008. — 256 с., ил.

Учебное пособие предназначено для практического освоения языка
FreeFem++, используемого для численного решения задач математической
физики методом конечных элементов. Содержит подробные программные
коды и большой набор задач. Рекомендуется студентам факультета математики, механики и компьютерных наук ЮФУ.

c⃝
М. Ю. Жуков, Е. В. Ширяева, 2008

c⃝
Издательство ЮФУ

Оглавление

Часть I. Основы метода конечных элементов
11

1
Метод конечных элементов в одномерном случае
12

1.1
Сильное и слабое решение задачи . . . . . . . . . . . . . . . . . . . . . . .
12

1.2
Построение приближенного решения задачи . . . . . . . . . . . . . . . . .
14

1.2.1
Слабое решение
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14

1.2.2
Метод Галеркина. Сильное решение . . . . . . . . . . . . . . . . . .
15

1.2.3
Конечно-разностный метод . . . . . . . . . . . . . . . . . . . . . . .
16

1.3
Выбор базисных функций
. . . . . . . . . . . . . . . . . . . . . . . . . . .
17

1.3.1
Финитные базисные функции
. . . . . . . . . . . . . . . . . . . . .
19

1.3.2
Вычисление элементов матрицы Aik и вектора bi . . . . . . . . . .
22

1.4
Естественные и главные краевые условия
. . . . . . . . . . . . . . . . . .
24

2
Метод конечных элементов в двумерном случае
29

2.1
Задача Дирихле для уравнения Лапласа. Слабое решение . . . . . . . . .
30

2.2
Построение приближенного решения . . . . . . . . . . . . . . . . . . . . .
31

Часть II. Обучение на примерах
35

3
Решение задач для уравнения Лапласа
38

3.1
Задача о стационарном распределении температуры . . . . . . . . . . . .
38

3.1.1
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38

3.1.2
Слабая формулировка задачи . . . . . . . . . . . . . . . . . . . . .
39

3.1.3
Слабая формулировка задачи на языке FreeFem++ . . . . . . . . .
41

3.1.4
Задание области D на языке FreeFem++ . . . . . . . . . . . . . . .
42

3.1.5
Полный код программы на языке FreeFem++ . . . . . . . . . . . .
43

3.2
Решение задачи о распределении температуры в областях сложной формы 44
3.2.1
Решение задачи для круга . . . . . . . . . . . . . . . . . . . . . . .
45

3.2.2
Распределение температуры в четырехугольнике . . . . . . . . . .
45

3.2.3
Решение задачи в криволинейной области . . . . . . . . . . . . . .
46

3.2.4
Решение задачи в области с отверстием
. . . . . . . . . . . . . . .
47

3.3
Физические задачи, приводящие к уравнению Лапласа . . . . . . . . . . .
48

3.3.1
Теплопроводность . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48

3.3.2
Диффузия
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51

3.3.3
Электрический потенциал (электростатика) . . . . . . . . . . . . .
52

3.3.4
Электрический потенциал (проводимость) . . . . . . . . . . . . . .
53

3.3.5
Потенциальное течение несжимаемой жидкости . . . . . . . . . . .
54

3.3.5.1
Стационарное обтекание крыла . . . . . . . . . . . . . . .
55

4
Решение нестационарной задачи для уравнения Лапласа
59

4.1
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59

4.2
Способ построения решения
. . . . . . . . . . . . . . . . . . . . . . . . . .
60

4.2.1
Аппроксимация производной по времени . . . . . . . . . . . . . . .
60

4.2.2
Слабая формулировка задачи . . . . . . . . . . . . . . . . . . . . .
61

3

Оглавление

4.3
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . . . . .
61

4.4
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . .
63

4.5
Общая схема аппроксимации производной по времени . . . . . . . . . . .
64

4.6
Контроль погрешности . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65

5
Уравнение переноса
69

5.1
Постановка задачи Коши для гиперболических уравнений . . . . . . . . .
69

5.2
Метод характеристик . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
70

5.3
Метод характеристик (двумерный случай) . . . . . . . . . . . . . . . . . .
72

5.4
Аппроксимация уравнения переноса . . . . . . . . . . . . . . . . . . . . . .
74

5.5
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . . . . .
75

5.6
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . .
76

5.7
Уравнение диффузии–переноса
. . . . . . . . . . . . . . . . . . . . . . . .
77

6
Уравнения реакция–диффузия. Окраска шкур животных
78

6.1
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78

6.2
Слабая формулировка задачи
. . . . . . . . . . . . . . . . . . . . . . . . .
79

6.3
Фрагменты кодов на языке FreeFem++ . . . . . . . . . . . . . . . . . . . .
80

6.4
Окраска шкур животных . . . . . . . . . . . . . . . . . . . . . . . . . . . .
81

6.4.1
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . .
82

6.4.1.1
Деформация области . . . . . . . . . . . . . . . . . . . . .
84

6.4.1.2
Различные модификации кода . . . . . . . . . . . . . . . .
84

6.4.2
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . .
85

7
Перенос-диффузия вихря
89

7.1
Переменные вихрь-функция тока
. . . . . . . . . . . . . . . . . . . . . . .
89

7.2
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90

7.3
Алгоритм решения задачи
. . . . . . . . . . . . . . . . . . . . . . . . . . .
91

7.4
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . . . . .
92

7.5
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . .
93

7.5.1
Движение вихрей в квадратной области
. . . . . . . . . . . . . . .
93

7.5.2
Движение вихрей в круге . . . . . . . . . . . . . . . . . . . . . . . .
94

7.5.3
Периодические краевые условия. Вихри на торе . . . . . . . . . . .
96

8
Тепловая конвекция
98

8.1
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
98

8.2
Алгоритм решения . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
99

8.2.1
Аппроксимация по времени . . . . . . . . . . . . . . . . . . . . . . . 100

8.2.2
Метод проекций . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

8.2.3
Решение задачи Неймана . . . . . . . . . . . . . . . . . . . . . . . . 102

8.2.4
Общая схема решения задачи
. . . . . . . . . . . . . . . . . . . . . 103

8.3
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . . . . . 104

8.4
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.4.1
Теплоизолированная прямоугольная область . . . . . . . . . . . . . 106

8.4.2
Теплоизолированный круг . . . . . . . . . . . . . . . . . . . . . . . 107

8.4.3
Тепловая конвекция в «чайнике»
. . . . . . . . . . . . . . . . . . . 107

8.4.4
Тепловая конвекция в «электрочайнике» . . . . . . . . . . . . . . . 108

8.5
Числа подобия . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

9
Различные течения жидкости
110

9.1
Задача о течении жидкости в полости
. . . . . . . . . . . . . . . . . . . . 110

9.1.1
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

9.1.2
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . 111

9.1.3
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . 112

9.2
Течение в канале . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
9.2.1
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

9.2.2
Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . 115

9.2.3
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . 117
9.2.3.1
Течение в канале T-образной формы . . . . . . . . . . . . 117

Оглавление
5

9.2.3.2
Обтекание препятствия в канале
. . . . . . . . . . . . . . 118

9.3
Обтекание тел жидкостью. Дорожка Кармана . . . . . . . . . . . . . . . . 119
9.3.1
Обтекание прямоугольного тела вязкой жидкостью . . . . . . . . . 119

9.3.2
Обтекание круглого тела вязкой жидкостью . . . . . . . . . . . . . 119

9.3.3
Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . 120

10 Перенос пассивной примеси
122

10.1 Задача о движении пассивной примеси в жидкости . . . . . . . . . . . . . 122
10.2 Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
10.3 Алгоритм решения задачи на языке FreeFem++ . . . . . . . . . . . . . . . 123
10.4 Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 126

11 Перенос примеси электрическим полем
128

11.1 Основные уравнения . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
11.2 Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
11.3 Алгоритм решения . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
11.4 Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . . . . . 132
11.5 Бездиффузионная модель переноса . . . . . . . . . . . . . . . . . . . . . . 133

11.5.1 Условия на разрыве . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
11.5.2 Задача о распаде начального разрыва
. . . . . . . . . . . . . . . . 134

11.6 Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 136

11.6.1 Случай α > 0 (увеличение σ с ростом c) . . . . . . . . . . . . . . . 136
11.6.2 Влияние параметра α на перенос примеси . . . . . . . . . . . . . . 137

11.7 Метод решения задачи, не использующий оператор convect . . . . . . . . 138

12 Задача о движении двух примесей
140

12.1 Основные уравнения и постановка задачи
. . . . . . . . . . . . . . . . . . 140

12.2 Реализация алгоритма решения на языке FreeFem++ . . . . . . . . . . . . 141
12.3 Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 142
12.4 Решение задачи с использованием ключевого слова convect . . . . . . . . 143

13 Перенос примеси электрическим полем и жидкостью
145

13.1 Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
13.2 Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 146

13.2.1 Перенос примеси жидкостью и электрическим полем . . . . . . . . 146
13.2.2 Сравнение различных типов переноса . . . . . . . . . . . . . . . . . 148

14 Решение задачи со свободной границей
149

14.1 Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
14.2 Код на языке FreeFem++ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
14.3 Результаты расчетов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

15 Течение Куэтта-Тейлора между вращающимися цилиндрами
154

15.1 Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

15.1.1 Функция тока
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

15.2 Алгоритм решения . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

15.2.1 Аппроксимация по времени . . . . . . . . . . . . . . . . . . . . . . . 156
15.2.2 Проекционный метод . . . . . . . . . . . . . . . . . . . . . . . . . . 157
15.2.3 Слабая формулировка задачи . . . . . . . . . . . . . . . . . . . . . 157

15.3 Реализация алгоритма на языке FreeFem++ . . . . . . . . . . . . . . . . . 158
15.4 Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 160

Часть III. Конструкции языка FreeFem++ (краткий обзор)
161

Оглавление

16 Синтаксис
162

16.1 Типы данных . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162

16.1.1 Основные типы данных . . . . . . . . . . . . . . . . . . . . . . . . . 164
16.1.2 Глобальные переменные . . . . . . . . . . . . . . . . . . . . . . . . . 165

16.2 Системные команды . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
16.3 Математические операции
. . . . . . . . . . . . . . . . . . . . . . . . . . . 166

16.4 Функции одной переменной . . . . . . . . . . . . . . . . . . . . . . . . . . . 168

16.4.1 Встроенные функции . . . . . . . . . . . . . . . . . . . . . . . . . . 168

16.4.1.1 Элементарные функции
. . . . . . . . . . . . . . . . . . . 168

16.4.1.2 Случайные функции
. . . . . . . . . . . . . . . . . . . . . 170

16.4.1.3 Дополнительные математические функции
. . . . . . . . 170

16.5 Функции двух переменных . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

16.5.1 Задание функций с помощью формул . . . . . . . . . . . . . . . . . 170
16.5.2 Конечноэлементные функции
. . . . . . . . . . . . . . . . . . . . . 171

16.6 Оператор условного перехода
. . . . . . . . . . . . . . . . . . . . . . . . . 172

16.7 Циклы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172

16.7.1 Цикл с параметром . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
16.7.2 Цикл с условием . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

16.8 Операторы ввода/вывода. Файлы . . . . . . . . . . . . . . . . . . . . . . . 174
16.9 Массивы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175

16.9.1 Одномерные массивы . . . . . . . . . . . . . . . . . . . . . . . . . . 175
16.9.2 Двумерные массивы с целочисленными индексами . . . . . . . . . 178

16.10Разреженные матрицы
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180

16.10.1Создание разреженной матрицы . . . . . . . . . . . . . . . . . . . . 180
16.10.2Операции над разреженными матрицами . . . . . . . . . . . . . . . 186

17 Генерация сеток
191

17.1 Простейшая область. Ключевое слово square . . . . . . . . . . . . . . . . 191
17.2 Ключевое слово border . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192
17.3 Запись/чтение сгенерированных сеток
. . . . . . . . . . . . . . . . . . . . 194

17.4 Ключевое слово «triangulate» . . . . . . . . . . . . . . . . . . . . . . . . 195
17.5 Ключевое слово «movemesh» . . . . . . . . . . . . . . . . . . . . . . . . . . 196
17.6 Ключевое слово «adaptmesh» . . . . . . . . . . . . . . . . . . . . . . . . . . 197

18 Конечные элементы
201

18.1 Интерполяция кусочно-линейными функциями . . . . . . . . . . . . . . . 201

18.1.1 Барицентрические координаты . . . . . . . . . . . . . . . . . . . . . 202
18.1.2 Интерполяция на треугольнике
. . . . . . . . . . . . . . . . . . . . 204

18.1.3 Линии уровня
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205

18.1.4 FE-функции
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206

18.1.5 Кусочно-линейная интерполяция
. . . . . . . . . . . . . . . . . . . 207

18.2 Базисные функции в FreeFem++ . . . . . . . . . . . . . . . . . . . . . . . . 209

18.2.1 P0 элементы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
18.2.2 P1 элементы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
18.2.3 P2 элементы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
18.2.4 P1nc элементы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
18.2.5 Элементы P1b и P2b
. . . . . . . . . . . . . . . . . . . . . . . . . . 213

18.3 Векторнозначные FE-функции . . . . . . . . . . . . . . . . . . . . . . . . . 214

18.3.1 RT0 элементы
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215

18.4 Загружаемые конечные элементы . . . . . . . . . . . . . . . . . . . . . . . 217
18.5 Численное интегрирование . . . . . . . . . . . . . . . . . . . . . . . . . . . 218

18.5.1 Интегрирование по границе области
. . . . . . . . . . . . . . . . . 218

18.5.2 Интегрирование по области
. . . . . . . . . . . . . . . . . . . . . . 220

Оглавление
7

19 Способы записи и решения задач в языке FreeFem++
223

19.1 Ключевые слова problem и solve . . . . . . . . . . . . . . . . . . . . . . . . 224

19.1.1 Слабая форма задачи и краевые условия . . . . . . . . . . . . . . . 224
19.1.2 Параметры, влияющие на решение задачи . . . . . . . . . . . . . . 226

19.2 Вариационные формы и разреженные матрицы . . . . . . . . . . . . . . . 227
19.3 Ключевое слово macro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230

20 Задачи на собственные значения
233

20.1 Функция EigenValue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233
20.2 Собственные значения оператора Лапласа . . . . . . . . . . . . . . . . . . 234

20.2.1 Задача для прямоугольника . . . . . . . . . . . . . . . . . . . . . . 235
20.2.2 Задача для квадрата
. . . . . . . . . . . . . . . . . . . . . . . . . . 235

20.3 Вычислительный эксперимент . . . . . . . . . . . . . . . . . . . . . . . . . 236

21 Визуализация результатов расчетов
238

21.1 Визуализация с помощью средств FreeFem++ . . . . . . . . . . . . . . . . 238

21.1.1 Параметры команды plot
. . . . . . . . . . . . . . . . . . . . . . . 239

21.2 Визуализация с помощью программы medit
. . . . . . . . . . . . . . . . . 243

21.3 Визуализация с помощью программы gnuplot
. . . . . . . . . . . . . . . . 244

A Установка программного обеспечения
246

A.1 Взаимодействие с текстовыми редакторами
. . . . . . . . . . . . . . . . . 246

B Используемые обозначения и формулы
248

C Список ключевых слов
249

Литература
250

Предисловие

FreeFem++ это интегрированная среда разработки (IDE) со своим собственным высокоуровневым языком программирования, предназначенная
для численного решения дифференциальных уравнений в частных производных методом конечных элементов. FreeFem++ позволяет исследовать
стационарные и нестационарные математические модели (в физике, химии, биологии, инженерных приложениях и т. п.) в 2D-пространственном
случае и является прекрасным инструментом как для обучения методу
конечных элементов, так и для решения научных проблем.

FreeFem++ — свободно распространяемое математическое обеспечение,
которое используется на различных платформах (Unix OS, Windows 9x,
2000, NT, XP и MacOS X). Авторы FreeFem++ — F. Hecht, O. Pironneau,
A. Le Hyaric, K. Ohtsuka — сотрудники лаборатории Ж.-Л. Лионса университета Пьера и Марии Кюри (Париж VI) и Французской Академии Наук.
В комплект установки FreeFem++1, помимо самой интегрированной среды
разработки, включено большое количество примеров и подробное руководство пользователя [1]. Заметим, что хотя FreeFem++ и является средой
разработки, по существу, это все-таки компилятор специального языка,
существенно ориентированного на работу с методом конечных элементов
и поэтому в дальнейшем удобно говорить именно о языке FreeFem++.

Коротко перечислим возможности FreeFem++:
1. Построение численного решения стационарных и нестационарных,
линейных и нелинейных двумерных краевых задач.

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

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

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

5. Для триангуляции области используется встроенный автоматический
генератор сеток, основанный на алгоритме Делоне–Вороного, при этом
плотность внутренних точек области пропорциональна плотности точек

1 http://www.freefem.org/ff++

8

Предисловие
9

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

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

7. Для проведения расчетов имеется большой набор быстрых алгоритмов для решения эллиптических конечно-разностных задач. Это разнообразные прямые и итеративные алгоритмы решения систем линейных уравнений — LU, Cholesky, Crout, CG, GMRES, UMFPACK, а также алгоритмы
для нахождения собственных значений и собственных векторов.

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

9. Алгоритмы записываются на языке, близком к C++ (знание C++ не
требуется). При этом эффективность кода по скорости близка к оптимальной, сравнимой со скоростью программ, непосредственно написанных на
языке С++.

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

По сравнению с коммерческими пакетами, предназначенными для решения уравнений в частных производных, например, FemLab или FlexPDE, использование FreeFem++ требует б´ольшей «умственной» работы от пользователя. В частности, невозможно записать задачу в исходном виде (FlexPDE
предоставляет такую возможность, а FemLab содержит список задач, в котором можно задавать численные значения параметров) — требуется преобразование к так называемой слабой формулировке. Менее развит и интерфейс среды FreeFem++. Однако эти минусы почти полностью компенсируются простотой языка FreeFem++, возможностью доступа ко всем
внутренним данным и возможностью создания собственных алгоритмов.

Иными словами, FreeFem++ это не «черный ящик», решающий задачи,
а инструмент для исследовательской работы, которая, зачастую, требует конструирования уникальных методов решения. Одновременно с этим
простота языка FreeFem++, как уже говорилось, позволяет использовать
его для обучения способам применения метода конечных элементов.

Предлагаемая книга предназначена читателям, желающим использовать метод конечных элементов для решения своих задач. В ней содержится большое количество примеров применения средств языка FreeFem++ к
решению разнообразных задач математической физики. Книга ни в коем
случае не заменяет руководства по использованию FreeFem++ [1]. Более
того, предполагается, что чтение данной книги и [1] будет осуществляться
параллельно. В [1] читатель найдет большое количество примеров и сможет изучить практически все возможности языка FreeFem++.

Коротко опишем структуру материала в книге, который разбит на три

Предисловие

части. В первой части (гл. 1, 2) дано краткое изложение метода конечных
элементов для одномерного и двумерного случаев. Вторая часть — это
алгоритмы решения на языке FreeFem++ различных конкретных задач.
Спектр этих задач весьма широк: простейшие задачи для уравнения Лапласа в гл. 3, 4, задачи для уравнения переноса в гл. 5, задача об окраске
шкур животных в гл. 6. Именно в этих главах содержатся основные алгоритмы и приемы, которые затем активно используются. Гл. 7–15 содержат
задачи, связанные с течением жидкости: вихревые течения в гл. 7, тепловая конвекция в гл. 8, течения в каналах и обтекание препятствий в гл. 9,
течение жидкости со свободной границей в гл. 14, течение Куэтта-Тейлора
в гл. 15. Гл. 10–13 содержат задачи о переносе примесей в жидкости.

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

Третья часть содержит детальное описание языка FreeFem++ и предназначена для углубленного изучения его возможностей. В первую очередь
это описание синтаксиса в гл. 16, описание приемов управления триангуляцией области в гл. 17, описание принципов выбора базисных функций для
построения приближенного решения в гл. 18. В этой части, в гл. 19, описан
выбор методов численного решения эллиптических конечно-разностных
задач (систем линейных алгебраических уравнений) и изложены способы решения краевых задач, отличные от приведенных во второй части.
В гл. 20 изложен способ решения краевых задач на собственные значения
при помощи языка FreeFem++. Наконец, гл. 21 посвящена возможностям
визуализации результатов расчета.

Неоценимую помощь авторам оказали студенты факультета математики, механики и компьютерных наук Южного федерального университета Ю. К. Самадова, Д. Д. Дуброва и И. И. Чернухина (подготовка части
рисунков и тестирование алгоритмов), а также сотрудники университета
Н. М. Жукова, прочитавшая рукопись и лишившая будущего читателя возможности обнаружить большое количество опечаток, И. В. Ширяева, сделавшая ряд важных замечаний о синтаксисе языка при написании гл. 16.

Работа выполнена благодаря поддержке внутренних грантов Южного
федерального университета, а также при поддержке грантов РФФИ (07-0100389а), Европейского научного объединения «Регулярная и хаотическая
гидродинамика» (07-01-9213) и гранта Президента поддержки ведущих научных школ Российской Федерации (НШ.5747.2006.1).

Часть I
Основы метода конечных элементов

Метод конечных элементов (МКЭ) является популярным и эффективным методом численного решения различных задач математической физики. Привести список огромного количества литературы, посвященной
основам метода и проблемам, связанным с его применением, — практически невыполнимая задача. Метод давно стал классическим и в настоящее
время его описание традиционно включается во многие учебники по численным методам [2–5] (см. также [6–14]).

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

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

Скажем несколько слов об используемой терминологии. Несмотря на
то, что метод конечных элементов известен очень давно, до сих пор в литературе используются различные термины. Так, сам метод называется —
метод конечных элементов, вариационно-разностный метод, проекционносеточный метод и т. д. Для названия базисных функций используются термины: тестовые функции, пробные функции и т. п. Это связано с широкой
применимостью метода для решения различного рода задач — задач теории упругости, гидродинамики и переноса, электродинамики, инженерных
задач и пр. В каждой из упомянутых областей принята своя специфическая терминология, зачастую не очень математически строгая.

Глава 1

Метод конечных элементов
в одномерном случае

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

1.1
Сильное и слабое решение задачи

Рассмотрим первую краевую задачу (задачу Дирихле) для определения
функции u(x)

− d

dx

p(x)du(x)

dx

+ q(x)u(x) = f(x),
0 < x < 1,
(1.1)

u(0) = 0,
u(1) = 0,
(1.2)

где p(x), q(x), f(x) — известные функции.

Предположим, что функции q(x), f(x) непрерывны, а функция p(x) —
непрерывно дифференцируема. Решение задачи (1.1), (1.2), являющееся
дважды непрерывно дифференцируемой функцией, будем называть сильным решением.

Можно попытаться переформулировать задачу (1.1), (1.2) с целью снижения требований, накладываемых на непрерывность функций p(x), q(x),
f(x) и функции u(x). Для этого умножим уравнение (1.1) на некоторую

12

Глава 1. Метод конечных элементов в одномерном случае
13

функцию v(x) и проинтегрируем на отрезке [0, 1]

−

10

d
dx

p(x)du(x)

dx

v(x) dx +

10
q(x)u(x)v(x) dx =

10
f(x)v(x) dx.

Используя формулу интегрирования по частям, имеем

10
pdu

dx
dv
dx dx − pdu

dxv
1

0
+

10
quv dx =

10
fv dx.

Потребуем, чтобы функция v(x) удовлетворяла тем же краевым условиям, что и u(x), т. е.
v(0) = 0,
v(1) = 0.
(1.3)

Тогда

10

pdu

dx
dv
dx + quv − fv
dx = 0,
v(0) = 0,
v(1) = 0.
(1.4)

Если соотношение (1.4) выполнено для любых v(x), то оно называется
слабой формой (или вариационной формой) записи задачи (1.1), (1.2). На
самом деле, следует дополнительно указать какому классу принадлежат
функции v(x), например, v ∈ C1.

Сам способ получения соотношения (1.4) показывает, что если функция u(x) является решением задачи (1.1), (1.2), то эта же функция будет
решением (1.4). Обратное же утверждение в общем случае неверно. Это
видно уже из того, что для функций p(x), q(x), f(x) и u(x), входящих в
(1.4), можно требовать гораздо меньших ограничений на гладкость, чем
для соответствующих функций задачи (1.1), (1.2). Так, например, можно
считать функции p(x), q(x), f(x) непрерывными (p, q, f ∈ C0), а функции
u(x), v(x) непрерывно дифференцируемыми (u, v ∈ C1). На самом деле,
можно даже требовать от p(x), q(x), f(x) кусочной непрерывности, а от
u(x), v(x) — кусочной гладкости. Заметим также, что для функции u(x),
входящей в (1.4), не требуется даже выполнения краевых условий (1.2).

Функцию u(x), являющуюся решением задачи (1.4), будем называть
слабым решением задачи (1.1), (1.2). Из вышесказанного ясно, что сильное и слабое решения задачи (1.1), (1.2) в общем случае не совпадают.
Подчеркнем, что сильные и слабые решения могут различаться и в случае, когда для функций p(x), q(x), f(x) сохранены такие же требования
на гладкость, как в исходной задаче (1.1), (1.2).

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

1 В частности, FreeFem++ предназначен для решения именно таких задач.

1.2. Построение приближенного решения задачи

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

Далее вопросы, связанные со сходимостью слабых решений к сильному,
не рассматриваются, т. к. основная цель книги — это описание работы с
FreeFem++. Ограничимся лишь ссылкой на [3, 5], где эти вопросы исследуются применительно к методу конечных элементов и в которых имеется
достаточно обширная библиография.

1.2
Построение приближенного решения задачи

1.2.1
Слабое решение

Будем искать приближенное решение задачи (1.4) в виде

uh(x) =

n−1
k=1
ckϕk(x),
(1.5)

где uh(x) — приближенное решение, функции ϕk(x) предполагаются известными, линейно независимыми и называются базисными функциями
(другие названия — тестовые функции, пробные функции), ck — коэффициенты, подлежащие определению.

Формула (1.5) задает аппроксимацию функции u(x) в виде некоторого сужения (проектирования) на конечномерное пространство, определяемое базисом ϕk(x), k = 1, . . . , n − 1. Точность аппроксимации, естественно, будет зависеть от того, насколько «хорошо» выбранное конечномерное
пространство приближает исходное пространство, которому должна принадлежать функция u(x). Заметим, что величина uh при заданных ϕk(x)
полностью определяется набором чисел (c1, c2, . . . , cn−1), т. е. можно считать, что uh = (c1, c2, . . . , cn−1) является вектором, принадлежащим Rn−1.

Подставляя (1.5) в (1.4), получим

10

p

n−1
k=1
ckϕ′
kv′ + q

n−1
k=1
ckϕkv − fv

dx = 0.

Учитывая, что (1.4) должно выполняться для любых v(x), принадлежащих некоторому классу, выберем в качестве v(x) некоторый набор функций, а именно v(x) = ϕi(x), i = 1, . . . , n−1. В соответствии с требованиями
v(0) = 0, v(1) = 0 (см. (1.4)) будем считать, что

ϕk(0) = 0,
ϕk(1) = 0,
k = 1, . . . , n − 1.
(1.6)

В этом случае имеем

10

p

n−1
k=1
ckϕ′
kϕ′
i + q

n−1
k=1
ckϕkϕi − fϕi

dx = 0,
i = 1, . . . , n − 1.

Глава 1. Метод конечных элементов в одномерном случае
15

Перепишем это соотношение в виде

n−1
k=1
Aikck = bi,
i = 1, . . . , n − 1,
(1.7)

где введены следующие обозначения

Aik =

10

pϕ′
iϕ′
k + qϕiϕk
dx,
bi =

10
fϕi dx,
i, k = 1, . . . , n − 1.
(1.8)

Таким образом, для определения ck имеем систему линейных алгебраических уравнений. Решая эту систему и подставляя ck в (1.5), получим
приближенное решение задачи (1.4). Подчеркнем, что приближенное решение (1.5) в силу условий (1.6) автоматически удовлетворяет краевым
условиям (1.2) для исходной задачи.

1.2.2
Метод Галеркина. Сильное решение

Разыскивать решение в виде (1.5) можно и для исходной задачи. Подставляя (1.5) в (1.1), имеем

δ(x) = − d

dx

p

n−1
k=1
ckϕ′
k

+ q

n−1
k=1
ckϕk − f

или

δ(x) = −p′
n−1
k=1
ckϕ′
k − p

n−1
k=1
ckϕ′′
k + q

n−1
k=1
ckϕk − f.
(1.9)

Здесь δ(x) — невязка, возникающая после подстановки приближенного решения uh(x) в уравнение (1.1).

Умножим соотношение (1.9) на ϕi и, проинтегрировав на [0, 1], потребуем выполнения равенств

10
δ(x)ϕi(x) dx = −

10

n−1
k=1
ck (p′ϕ′
k + pϕ′′
k − qϕk) + f

ϕi dx = 0
(1.10)

или

n−1
k=1
ck

10
(−p′ϕ′
kϕi − pϕ′′
kϕi + qϕkϕi) dx −

10
fϕi dx = 0,
i = 1, . . . , n − 1.

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

n−1
k=1
˜Aikck = bi,
i = 1, . . . , n − 1,