Численное решение краевых задач

Министерство науки и высшего образования Российской Федерации 

НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ 

 
 
 
 
 
 
Ю.А. ГОСТЕЕВ 
 
 
 
 
 
 
ЧИСЛЕННОЕ РЕШЕНИЕ  
КРАЕВЫХ ЗАДАЧ 
 
Учебно-методическое пособие 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
НОВОСИБИРСК 
2019 

 

УДК 53:519.6(075.8) 
          Г 725 
 
 

Рецензенты: 

д-р техн. наук, профессор А.В. Чичиндаев 
канд. техн. наук, доцент Ю.В. Телкова 
 
Работа подготовлена на кафедре аэрогидродинамики 
для студентов III и IV курсов ФЛА направления «Баллистика  
и гидроаэродинамика» и утверждена Редакционно-издательским  
советом университета в качестве учебно-методического пособия 
 
 
Гостеев Ю.А. 
Г 725   
Численное решение краевых задач: учебно-методическое 
пособие / Ю.А. Гостеев. – Новосибирск: Изд-во НГТУ, 2019. – 
66 с. 

ISBN 978-5-7782-4075-9 

На примерах рассмотрены основные конечно-разностные методы 
решения краевых и начально-краевых задач для уравнений математической физики. Приведены расчетные программы на языке Фортран и 
задания расчетно-графических работ. Показаны приемы работы в пакете прикладных программ при решении краевых задач методом конечных объемов. В приложении содержится необходимый материал 
для создания графических пользовательских приложений. 
Пособие предназначено для студентов III и IV курсов ФЛА направления «Баллистика и гидроаэродинамика», изучающих дисциплины 
«Вычислительная математика» и «Численные методы механики жидкости и газа». 
 
 
УДК 53:519.6(075.8) 
 
 
ISBN 978-5-7782-4075-9  
 
 
 
 
© Гостеев Ю.А., 2019 
© Новосибирский государственный  
    технический университет, 2019 

 

ВВЕДЕНИЕ 

В пособии на примерах рассмотрены основные конечно-разностные 
методы решения краевых и начально-краевых задач для уравнений математической физики. Приведены расчетные программы на языке Фортран и задания расчетно-графических работ. Показаны приемы работы в 
пакете прикладных программ при решении краевых задач методом конечных объемов. В приложении содержится необходимый материал для 
создания графических пользовательских приложений. 
Успешное усвоение материала пособия предполагает знание основ 
математического анализа, специальных разделов высшей математики, 
программирования на алгоритмическом языке Фортран. 
Пособие предназначено для студентов III и IV курсов ФЛА направления «Баллистика и гидроаэродинамика», изучающих дисциплины 
«Вычислительная математика» и «Численные методы механики жидкости и газа». 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 

1. ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ  
ДЛЯ ОДУ  

ЦЕЛЬ РАБОТЫ 

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

ЗАДАНИЕ 

Методами конечных разностей (МКР) и конечных объемов (МКО) 
решить краевую задачу для поля температуры плоского слоя: 

2

2
( )
0,
0
d T
S x
x
l
dx




 , 

 
л
T
T

   при   
0
x 
, 
(1) 

п
T
T

   при   x
l
 . 

Здесь x – координата, l  – толщина слоя;   – коэффициент теплопроводности материала; S  – источник тепла; 
л
T  и 
п
T  – граничные температуры. Определить потоки тепла на границах 
л,п
0,
(
/
) x
L
q
dT dn

 
, 

где n – внешняя нормаль. 
Варианты заданий приведены в таблице. 

Номер 
варианта Материал 
l , м 
л
T ,С 
п
T ,С 
S , МВт/м3 

1 
Медь 
0,05 
200 
10 
cos(
/ )
x l

 

2 
Сталь 
0,10 
20 
300 
2
/
x l  

О к о н ч а н и е  т а б л и ц ы  

Номер 
варианта Материал 
l , м 
л
T ,С 
п
T ,С 
S , МВт/м3 

3 
Титан 
0,15 
500 
10 
3
2 /x l
e
 

4 
Никель 
0,20 
100 
600 
4 cos(5
/ )
x l

 

5 
Медь 
0,25 
50 
100 
5 3
/
x l  

6 
Сталь 
0,20 
150 
10 
4
/x l
e
 

7 
Титан 
0,15 
15 
20 
3 sin (2
/ )
x l

 

8 
Никель 
0,10 
400 
15 
2 4
/
x l  

9 
Медь 
0,05 
30 
350 
3 /x l
e
 

10 
Сталь 
0,10 
200 
20 
2 cos(3
/ )
x l

 

11 
Титан 
0,15 
10 
250 
3 5
/
x l  

12 
Никель 
0,20 
300 
30 
4 sin (0,5
/ )
x l

 

13 
Медь 
0,24 
450 
20 
3 1 3 /
x l

 

14 
Сталь 
0,22 
120 
570 
4 1 3 /
x l
e 
 

15 
Титан 
0,16 
60 
90 
5 cos(1
4
/ )
x l
 
 

16 
Никель 
0,11 
140 
20 
3 3 2
/
x l

 

17 
Медь 
0,07 
30 
230 
6 1 2 /x l
e 
 

18 
Сталь 
0,12 
330 
25 
2 sin (1
2
/ )
x l
 
 

19 
Титан 
0,14 
50 
340 
3 4 1
/
x l

 

20 
Никель 
0,18 
220 
15 
2 1
/x l
e 
 

УКАЗАНИЯ 

Точное решение задачи (1) имеет вид 



л
2
п
л
2
2
2
( )
(0)
(0)
( )
( )
e
x
T x
T
r
T
T
r
r l
r x
l







, 

 
2( )
( )
r
x
dx r x dx
 

,   ( )
( ) /
r x
S x
 
 . 
(2) 

Дифференцируя (2), найдем точные выражения для 
л,п
q
. 

МКР 

Задача (1) есть частный случай линейной краевой задачи для ОДУ 
второго порядка 

 
л
п
( )
( )
( ),
,
( )
,
( )
u
p x u
q x u
r x
a
x
b
u a
u
u b
u









. 
(3) 

1. На отрезке [ , ]
a b  введите разностную сетку 
1
2
(
,
, ...,
)
N
x x
x
, где 
(
1)
ix
a
i
h



, 
1,...,
,
(
) / (
1)
i
N h
b
a
N




. С помощью формул численного дифференцирования замените дифференциальную задачу (3) 
разностной: 

1
л
u
u

, 

 
1
1
1
1
2
2
(
)
( )
(
),
2,...,
1
2

i
i
i
i
i
i
i
i
u
u
u
u
u
p x
q x u
r x
i
N
h
h












 , 
(4) 

п
N
u
u

. 

Собирая коэффициенты при 
1
iu  , 
iu , 
1
iu  , найдите для вектора неизвестных ( 1u ,…,
N
u
) СЛАУ с трехдиагональной матрицей, которую решаем методом прогонки. Для определения граничных потоков тепла 1q   
и 
N
q   используйте аппроксимацию производной по несимметричным 
формулам второго порядка: 

1
2
3

1

3
4
2
u
u
u
du
dx
h




 




,     
2
1
4
3
2

N
N
N

N

u
u
u
du
dx
h












. 

2. Напишите фортран-программу, реализующую метод конечных 
разностей (4). За основу возьмите программу Kraev1 (см. ниже), которая 
решает задачу (1) для слоя Al (
1
1
202,4 Вт м К


 

) толщиной l  1 м с 
граничными температурами 
л
u  300 К, 
п
u  500 К и источником тепла 

max
( )
sin(5
)
S x
S
x


, 
max
S
 1,6192 МВт/м3. Данная задача имеет точное 
решение: 

max
max
л
п
л
2
2
( )
sin(5
)
sin(5
)
(5 )
(5 )
e
r
r
x
u
x
u
u
u
l
x
l
















,   max
max /
r
S

 . 

3. Используя программу, проведите расчеты для числа шагов 
(
1)
N 
 10, 20, 40 и 80. Полученные данные занесите в таблицу. 

N 
max

  maxue(xi) – ui 
qл – q1 
qп – qN 

 
 
 
 

На основании анализа погрешности сделайте вывод о порядке аппроксимации метода. Постройте графики точного и приближенного решения для разных N, убедитесь в сеточной сходимости приближенного 
решения (т. е. убедитесь в уменьшении влияния сетки на решение при 
ее измельчении). 
 
program kraev1 
parameter(N=11) 
dimension c(N),d(N),e(N),f(N),x(N),u(N),ue(N) 
common /par1/ xR,uL,uR,rmax 
open(10,file='res.dat') 
xL=0.0; uL=300 
xR=1.0; uR=500 
rmax=8e3 ! Smax/lambda 
h=(xR-xL)/(N-1); h2=h*h 
do j=1,N;  x(j)=xL+(j-1)*h; 
c(j)=1.0-(h/2)*p(x(j)); d(j)=h2*q(x(j))-2.0 
e(j)=1.0+(h/2)*p(x(j)); f(j)=h2*r(x(j)) 
end do 
c(1)=0.0; e(1)=0.0; d(1)=1.0; f(1)=uL 
c(N)=0.0; e(N)=0.0; d(N)=1.0; f(N)=uR 
call Progon(N,c,d,e,f,u) 
print*,'      x           u           uex' 
do j=1,N 
 ue(j)=uex(x(j)) 
 print '(3(1x,1pe12.5))',x(j),u(j),ue(j) 
 write(10,100)x(j),u(j),ue(j) 
end do 
du=maxval(abs(u-ue))/uL 
print*,' du=',du 
100 format(3(1x,1pe12.5)) 
end 

function p(x) 
p=0.0 
end 
function q(x) 
q=0.0 
end 
 
function r(x) 
common /par1/ xR,uL,uR,rmax 
pi=4.*atan(1.) 
r=-rmax*sin(5.*pi*x) 
end 
 
function uex(x) 
common /par1/ xR,uL,uR,rmax 
pi=4.*atan(1.); rmax5=rmax/(5*pi)**2 
uex=uL+((uR-uL)-rmax5*sin(5*pi*xR))*x/xR+rmax5*sin(5*pi*x) 
end 
 
Результаты: 
       x           u           uex 
  0.00000E+00  3.00000E+02  3.00000E+02 
  1.00000E-01  3.60000E+02  3.52423E+02 
  2.00000E-01  3.40000E+02  3.40000E+02 
  3.00000E-01  3.20000E+02  3.27577E+02 
  4.00000E-01  3.80000E+02  3.80000E+02 
  5.00000E-01  4.40000E+02  4.32423E+02 
  6.00000E-01  4.20000E+02  4.20000E+02 
  7.00000E-01  4.00000E+02  4.07577E+02 
  8.00000E-01  4.60000E+02  4.60000E+02 
  9.00000E-01  5.20000E+02  5.12423E+02 
  1.00000E+00  5.00000E+02  5.00000E+02 
  du=  2.5257669E-02 

МКО (ПАКЕТ ANSYS) 

1. Создание проекта (программа Workbench). После запуска Workbench раскройте в «Toolbox» закладку «Analyses Systems» и выберите 
«Fluid Flow (Fluent)» (двойной клик «мыши» либо перетаскивание  

на поле «Project Schematic»). Сохраните проект (кнопка «Save») в специально созданной рабочей папке. Путь к рабочей папке не должен содержать букв кириллицы. 

 
2. Создание геометрической модели расчетной области (DesignModeler). Выполните двойной клик «мыши» по пункту «Geometry» в 
схеме проекта. При первом запуске выберите единицу измерения 
длины («Meter» по умолчанию). Далее создайте на плоскости xy  чертеж: выделите «XYPlane», нажмите кнопку «New Sketch» 
; на закладке «Sketching» с помощью инструментов «Draw» ->  «Rectangle» 
и «Dimensions» -> «Horizontal», «Vertical» начертите прямоугольную 
область длиной l  ( 0
x
l

 ) и высотой / 40
l
 ( 0
/ 40
y
l


). (Для тестовой задачи l  1 м.) Затем переключитесь на закладку «Mode- 
ling», выберите пункт меню «Concept»->«Surfaces From Sketches» 
 

 

и нажмите кнопку 
. Сохраните проект и выйдите из DesignModeler. 
3. Генерация расчетной сетки и задание границ (Meshing). Выполните двойной клик «мыши» по пункту «Mesh» в схеме проекта. В древовидном списке выберите «Mesh» (ПКМ) -> «Insert» -> «Sizing». Выделите нижнюю и верхнюю горизонтальные границы области (при этом 
соответствующий фильтр выделения объектов должен быть активным: 

) и в окне свойств нажмите «Apply». Также в окне 
свойств задайте «Type» = «Number of Divisions», «Number of Divisions» = 10, «Behavior» = «Hard». Аналогичные действия выполните для 
левой и правой вертикальных границ, причем задайте число ячеек 
«Number of Divisions» = 1. Далее выполните «Mesh» -> «Insert» -> 
«Mapped Face Meshing» и в свойствах укажите область для «Geometry». 
Затем нажмите 
. Будет сгенерирована конечно-объемная сетка с ячейками прямоугольной формы. Теперь создайте граничные поверхности: выделяя граничную линию, вызывайте контекстное 
меню (ПКМ) -> «Create Named Selection». Левая граница – «left», правая – «right», нижняя – «bot», верхняя – «top». Перед сохранением про
екта выполните 
. 
4. Задание физической модели, начальных и граничных условий, параметров численного алгоритма (Fluent). Запуск препроцессора осуществляется с помощью кнопки «Setup» в дереве проекта. 
В пункте «General» раздела «Solution Setup» выберите опцию 
«Time» = «Steady» для решения стационарной задачи. 
В пункте «Models» подключите уравнение энергии «Energy – On». 
В пункте «Materials» выберите материал твердого тела: выделите 
«Solid» и нажмите кнопку «Create/Edit». В окне «Create/Edit Materials» 
вызовите базу данных «Fluent Database…», измените «Material Type» на 
«Solid» и выберите в списке требуемый материал («aluminum (al)» для 
тестовой задачи), затем нажмите «Copy» и «Close». В окне «Create/Edit 
Materials» выберите «Change/Create» и «Close». 
Создайте в каталоге проекта файл «source.с», содержащий пользовательскую функцию (UDF) для расчета источника тепла. Содержимое 
файла для тестовой задачи: