Применение численных методов для построения разностных моделей

МИНИСТЕРСТВО НАУКИ И ВЫСШЕГО ОБРАЗОВАНИЯ

РОССИЙСКОЙ ФЕДЕРАЦИИ

Федеральное государственное автономное 

образовательное учреждение высшего образования

«ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ»

Инженерно-технологическая академия

В. В. СЕМЕНИСТЫЙ

И. Э. ГАМОЛИНА
В. В. ДУРЯГИНА

ПРИМЕНЕНИЕ ЧИСЛЕННЫХ МЕТОДОВ 

ДЛЯ ПОСТРОЕНИЯ РАЗНОСТНЫХ МОДЕЛЕЙ

Учебное пособие

Ростов-на-Дону – Таганрог

Издательство Южного федерального университета

2021

Оглавление

2

УДК 518.12(075.8)
ББК 22.162 Я 73

С301

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

Института компьютерных технологий и информационной безопасности 

Южного федерального университета 
(протокол № 10 от 11 февраля 2020 г.)

Рецензенты:

кандидат педагогических наук,

заведующий кафедры теории и методики математического образования 

Института математики, механики и компьютерных наук 

им. И. И. Воровича Южного федерального университета Ю. В. Романов

доктор физико-математических наук, профессор,

заведующий отделом стихийных явлений 

ФГБУ «Высокогорный геофизический институт» А. Х. Аджиев

Семенистый, В. В.

С301
Применение численных методов для построения разностных моде-

лей : учебное пособие / В. В. Семенистый, И. Э. Гамолина, В. В. Дуря-
гина ; Южный федеральный университет. – Ростов-на-Дону ; Таганрог :
Издательство Южного федерального университета, 2021. – 119 с.

ISBN 978-5-9275-3765-5
Учебное пособие написано для магистрантов первого курса обучения и 

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

УДК 518.12(075.8)

ББК 22.162 Я 73

ISBN 978-5-9275-3765-5

© Южный федеральный университет, 2021
© Семенистый В. В., Гамолина И. Э.,

Дурягина В. В., 2021

© Оформление. Макет. Издательство

Южного федерального университета, 2021

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ ……………………………………………………………
5

1.  РАЗНОСТНЫЙ  МЕТОД  РАСЧЕТА  ЗАДАЧИ  СОУДАРЕНИЯ 
ГИБКИХ  ПЛАСТИН ………………………………………………….
7

1.1. Постановка задачи ……………………………………………..
7

1.2. Постановка дифференциальной задачи ………………………
9

1.3. Построение разностной схемы ………………………………..
12

1.4. Расчет параметров процесса …………………………………..
21

Вопросы для самоконтроля ………………………………………...
31

2.  РАЗНОСТНЫЙ  МЕТОД  РАСЧЕТА  ЗАДАЧИ  ЭЛЕКТРИЧЕ-
СКОГО  СОСТОЯНИЯ  ГОРИЗОНТАЛЬНО-ОДНОРОДНОГО  ТУР-
БУЛЕНТНОГО  ПРИЗЕМНОГО СЛОЯ ……………………………….
33

2.1. Постановка дифференциальной задачи ………………………
33

2.2. Построение сеточной модели …………………………………
35

2.2.1. Алгоритм исследования дифференциальной задачи для 
построения сеточной модели …………………………………………
35

2.2.2. Сеточное уравнение для первой модели ……………………..
36

2.2.3. Сеточное уравнение для второй модели ……………………..
38

2.2.4. Сеточное уравнение для третьей модели …………………...
38

2.2.5. Сеточное уравнение для четвертой модели ………………..
39

2.2.6. Решение уравнения напряженности электрического поля
40

2.2.7. Разностная схема для вычисления функции концентрации 
ионов ……………………………………………………………………….
41

2.2.8. Разностная модель исследуемой задачи в матричной форме
41

2.3. Основные положения теории и примеры решения задач …...
43

Вопросы для самоконтроля ………………………………………...
54

Варианты контрольных работ ……………………………………...
58

3. КРАЕВЫЕ  ЗАДАЧИ  ДЛЯ  ЭЛЛИПТИЧЕСКИХ  УРАВНЕНИЙ
62

3.1. Теоретическая часть ……………………………………………
62

3.2. Практическая часть …………………………………………….
63

4.  НЕСТАЦИОНАРНЫЕ ЗАДАЧИ МАТЕМАТИЧЕСКОЙ ФИ-
ЗИКИ …………………………………………………………………...
70

4.1. Теоретическая часть ……………………………………………
70

Оглавление

4

4.2. Практическая часть …………………………………………….
72

5.  КРАЕВЫЕ ЗАДАЧИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕН-
ЦИАЛЬНЫХ УРАВНЕНИЙ ………………………………………...
78

5.1. Теоретическая часть ……………………………………………
78

5.2. Практическая часть …………………………………………….
83

Вопросы для самоконтроля ………………………………………...
88

6. ЗАДАЧИ ДЛЯ САМОСТОЯТЕЛЬНОЙ РАБОТЫ …………….
91

6.1. Задачи для уравнений эллиптического типа …………………
91

6.2. Нестационарные краевые задачи ……………………………...
94

6.3. Задачи для обыкновенных дифференциальных уравнений …
97

Вопросы для самоконтроля…………………………………………
102

Проектное задание ………………………………………………….
105

ЗАКЛЮЧЕНИЕ ……………………………………………………….
106

СПИСОК  ЛИТЕРАТУРЫ ……………………………………………
107

ПРИЛОЖЕНИЯ. Варианты контрольных работ ……………………
108

ВВЕДЕНИЕ

В данном пособии излагаются основы численных методов, исполь-

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

В первой части работы (первые два раздела) уделяется внимание по-

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

В третьем, четвертом и пятом разделах учебного пособия прово-

дится исследование различных свойств разностных моделей.

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

применении различных граничных условий.

Если систематизировать стационарные задачи математической фи-

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

В данном учебном пособии представлены некоторые методы, в част-

ности итерационные, для решения сеточных уравнений.

Кроме того, рассмотрены разностные методы для получения числен-

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

Кроме того, рассмотрен вопрос, связанный с исследованием схем, 

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

Введение

6

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

В каждом разделе используется своя нумерация формул.
Данное учебное пособие может быть использовано для подготовки 

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

1.  РАЗНОСТНЫЙ МЕТОД РАСЧЕТА ЗАДАЧИ 

СОУДАРЕНИЯ ГИБКИХ ПЛАСТИН

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

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

удара – задачи о соударении двух сжимаемых пластин, уравнения которой 
описывается моделью идеальной жидкости [1].

В качестве физической модели пластины рассматривается тело та-

кой формы, у которой один из габаритных размеров, назовем его толщиной,
много меньше двух других габаритных размеров. В качестве основных ха-
рактеристик для пластин выберем следующие: а) срединную поверхность –
плоскость, равноудаленную от наружных поверхностей пластины; б) тол-
щину h; в) величину прогиба w, когда на пластину действуют внешние 
силы (нагрузка). Если провести классификацию по величине прогиба, то 
можно выделить следующие классы пластин. Первый вид пластин – жест-
кие пластины, в данном классе пластин величина прогиба w не может быть 
выше 20–25% от толщины h пластины. Если мы рассмотрим вид функцио-
нальной зависимости между прогибом и внешней нагрузкой в данном типе 
пластин, то выделим преобладание линейной зависимости, т.е. в основном 
возникают изгибные напряжения. Второй класс пластин – гибкие пла-
стины, в срединной плоскости которых возникают цепные или мембранные 
напряжения. Третий класс – это абсолютно гибкие пластины. В таком 
классе пластин функциональная зависимость между прогибом и нагрузкой 
является нелинейной. Доминируют только мембранные напряжения, т.е.
можно пренебречь изгибными напряжениями. 

По сравнению с одномерной нестационарной задачей о движении 

газа в трубе под действием поршня задача о соударении двух сжимаемых 
пластин более сложная. Усложнение связано, прежде всего, с появлением 
второй взаимодействующей деформируемой среды, а также с необходимо-
стью реализации более сложных граничных условий [2].

Подробно рассмотрим весь процесс постановки и решения задачи, 

а именно: постановку дифференциальной задачи, построение разностной 

1. Разностный метод расчета задачи соударения гибких пластин

8

модели, составление программы расчета параметров процесса взаимодей-
ствия пластин, результаты расчета и анализ особенностей процесса. При 
построении разностной модели используем метод Мейдера.

Рис.1. Расчетная область

Будем считать, что пластина 1 толщиной 𝛿1 из материала с началь-

ной плотностью 𝜌1 соударяется со скоростью 𝑢0 с неподвижной пласти-
ной 2, имеющей толщину 𝛿2 и начальную плотность 𝜌2 (рис. 1).

Неподвижную систему отсчета наблюдателя (эйлерова система от-

счета) зададим следующим образом. Выберем в качестве точки отсчета 
точку 0 в плоскости соударения. Движение обеих взаимодействующих пла-
стин обладает очевидной симметрией слоя: все текущие значения парамет-
ров движения и состояния среды (массовой скорости и, давления Р, плот-
ности 𝜌 и удельной внутренней энергии Е) будут одинаковыми в любой 
плоскости, параллельной плоскости соударения. Поэтому в качестве системы 
координат следует выбрать декартову прямоугольную систему (х, у, z), ось 
х которой ориентирована в направлении удара, а оси у и z лежат в плоскости 
соударения. При таком выборе системы координат вектор скорости имеет 
лишь одну компоненту 𝑣𝑥 = 𝑢, а параметры движения и состояния зависят 
лишь от одной координаты х и времени t. Это одномерная плоская неста-
ционарная задача.

1.2. Постановка дифференциальной задачи

9

1.2. Постановка дифференциальной задачи

Решение такой задачи целесообразно искать, описывая движение с 

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

За лагранжевы линейные координаты X индивидуальных точек 

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

Х = 𝑥|𝑡 = 0.
0
t
X
x =
=

Под лагранжевой массовой координатой т будем понимать погон-

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

𝑑𝑚 = 𝜌0𝑑𝑋,

где 𝑑𝑋 – соответствующий дифференциал лагранжевой линейной координаты (
начальное расстояние между двумя бесконечно близкими индивидуальными 
точками, между которыми заключена данная индивидуальная частица). 
При таком определении лагранжевы массовые координаты поверхностей 
пластин равны: на наружной поверхности пластины 1 (левая на 
рис. 1) 𝑚 = 0; на поверхности соударения (правая для пластины 1 и левая 
для пластины 2) 𝑚 = 𝜌01𝛿1; на наружной (правой) поверхности пластины 
𝑚 = 𝜌01𝛿1 + 𝜌02𝛿2. Будем использовать лагранжевы массовые координаты 
и искать решение в виде

и = и(т, t) 

Р = Р(т, t) 

𝜌 = 𝜌 (m, t) 

E = E(m, t) 

x = x(m, t), 

где х – изменяющиеся во времени эйлеровы координаты индивидуальных 
точек.

1. Разностный метод расчета задачи соударения гибких пластин

10

При использовании для материалов обеих пластин модели идеальной 

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

𝜕𝑢

𝜕𝑡 +

𝜕р

𝜕𝑚 = 0;
(1.1)

𝜕𝑥

𝜕𝑚 =

1

𝜌;
(1.2)

𝜕

𝜕𝑡 (𝐸 +

𝑢2

2 ) +

𝜕(𝑃𝑢)

𝜕𝑚 = 0;
(1.3)

𝜕𝑥

𝜕𝑡 = 𝑢;
(1.4)

𝑃 = 𝐴 [(

𝜌

𝜌0)

𝑛

− 1],
(1.5)

при этом константа в уравнении Тэта (1.5) зависит от свойств материала 
пластин и, следовательно, от координат:

𝐴 = 𝐴1, 𝜌0 = 𝜌01,  𝑛 = 𝑛1 при 0 ≤ 𝑚 ≤ 𝜌01𝛿1;
(1.6)

𝐴 = 𝐴2, 𝜌0 = 𝜌02,  𝑛 = 𝑛2 при 𝜌01𝛿1 ≤ 𝑚 ≤ 𝜌01𝛿1 + 𝜌02𝛿2.

Каждое дифференциальное уравнение практически явно выражает 

соответствующий закон сохранения. Очевидным своеобразием в этой си-
стеме обладают уравнение неразрывности (1.2) и уравнение энергии (1.3). 
Уравнение неразрывности записано в виде дифференциального соотноше-
ния между текущими эйлеровыми координатами х и лагранжевыми мас-
совыми координатами т индивидуальных точек среды, а уравнение энер-
гии – через удельную полную энергию Е + и2 /2.

Задача решается при следующих начальных условиях.
При 𝑡 = 0,  𝑢(𝑚, 0) = 𝑢0, 𝜌(𝑚, 0) = 𝜌01 для индивидуальных точек 

пластины 1 (0 ≤ 𝑚 ≤ 𝜌01𝛿1)  и u(m,0) = 0 , 𝜌 (m,0) = 𝜌02 для индивидуаль-
ных точек пластины 2 (𝜌01𝛿1 ≤ 𝑚 ≤ 𝜌01𝛿1 + 𝜌02𝛿2).

1.2. Постановка дифференциальной задачи

11

Начальные распределения давления и удельной внутренней энергии 

будем считать одинаковыми для обеих пластин:

Р = Р(m,0), E = Е(т,0)

при  0 ≤ 𝑚 ≤ 𝜌01𝛿1 + 𝜌02𝛿2.

Что касается начальных значений эйлеровых координат индивиду-

альных точек, то они совпадают со значениями лагранжевых линейных ко-
ординат:

х(т,0) = Х(т,0).

Граничные условия таковы. Будем считать, что на наружных свобод-

ных границах заданы поверхностные силы.

Для принятой модели идеальной жидкости динамические граничные 

условия на этих поверхностях сводятся к условиям:

на поверхностях обеих пластин т=0 и 0 ≤ 𝑚 = 𝑝01𝛿1 + 𝑝02𝛿2
отсутствует давление:

𝑝(0, 𝑡) = 0; 𝑝(𝑝01𝛿0 + 𝑝02𝛿2, 𝑡) = 0.
(1.7)

Граничные условия на обращенных друг к другу поверхностях пла-

стин (𝑚 = 𝜌01𝛿1) зависят от наличия или отсутствия контакта между ними 
(в общем случае возможны нарушение контакта в ходе взаимодействия 
пластин и переход к их свободному деформированию). О наличии или об 
отсутствии контакта можно судить по значениям эйлеровых координат ин-
дивидуальных точек обращенных друг к другу поверхностей пластин.

При 𝑥(𝜌01𝛿1 + 0) > 𝑥(𝜌01𝛿1 − 0) контакт отсутствует.
При 𝑥(𝜌01𝛿1 + 0) = 𝑥(𝜌01𝛿1 − 0) контакт имеется.
Случаю отсутствия контакта соответствуют динамические граничные 
условия на свободных поверхностях:

𝑝(𝜌01𝛿1 + 0) = 𝑝(𝜌01𝛿1 − 0) = 0,

при
𝑥(𝜌01𝛿1 + 0) > 𝑥(𝜌01𝛿1 − 0).
(1.8)

При наличии контакта между пластинами на поверхности контакта 

задаются граничные условия смешанного типа: условие «прилипания» по 
скоростям и условие равенства по давлениям, следующее из третьего закона 
Ньютона:

𝑢(𝜌01𝛿1 + 0) = 𝑢(𝜌01𝛿1 − 0), 𝑝(𝜌01𝛿1 + 0) = 𝑝(𝜌01𝛿1 − 0), 

1. Разностный метод расчета задачи соударения гибких пластин

12

при
𝑥(𝜌01𝛿1 + 0) = 𝑥(𝜌01𝛿1 − 0).
(1.9)

Таким образом, система уравнений (1.1)–(1.6) с приведенными выше 

начальными условиями и граничными условиями (1.7)–(1.9) математически 
описывает процесс соударения двух сжимаемых «жидких» пластин.

1.3. Построение разностной схемы

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

Псевдовязкость – это искусственный механизм, позволяющий осуществлять 
сквозной расчет ударных волн без явного выделения их фронтов 
на разностной сетке. Разностная схема с псевдовязкостью является однородной. 
Область фронта – область резкого изменения в распределениях параметров 
течения (при отсутствии сильных разрывов). Отметим, что следует 
различать ширину ударного перехода, обусловленную воздействием псевдо-
вязкости, и реальную ширину фронта ударной волны. В газах и жидкостях 
реальная ширина фронта ударной волны составляет несколько длин свободного 
пробега молекул. Эффективная же ширина «размазываемого» за 
счет псевдовязкости фронта определяется, по существу, параметрами сетки 
по координате. Псевдовязкость является аддитивной «квазивязкой» добав-
кой к давлению, т.е. сумма газодинамического давления и псевдовязкости опре-
деляет эффективное давление: 
=p+q. Поэтому модификация схемы с ис-

пользованием псевдовязкости сводится к изменению конечно-разностных 
уравнений движения (1.4) и энергии (1.3), в которых участвует давление 
(замена давления р эффективным давлением 
), а также непосредственно 

к определению сеточной функции псевдовязкости.

Принцип «фиктивной» ячейки состоит в том, что если слева от гра-

ницы ввести еще одну ячейку с массовым интервалом ℎ−1 = 0 и ассоциировать