Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
МЕТОДОЛОГИЯ НАУЧНОГО ТВОРЧЕСТВА
ЛЕКЦИЯ 13
ЛЕКЦИЯ 13
МЕТОДЫ РЕШЕНИЯ УЧП
Рассмотрим численное решение уравнений с частными производными (УЧП) на примере решения уравнения теплопроводности (диффузии)
где
- плотность тепловых источников, - коэффициент теплопроводности
Если , то мы имеем стационарное уравнение теплопроводности
- уравнение Пуассона.
Если плотность тепловых источников тоже равна 0, то получаемое уравнение называется уравнением Лапласа.
Стационарное уравнение теплопроводности может быть решено, если известны граничные условия.
Это могут быть или значения искомой функции на границе (задача Дирихле).
Или значения потока тепла на границе (задача Неймана)
Либо смешанные условия.
Если необходимо решать УЧП в области, имеющей круговую симметрию, оператор Лапласа удобна записать в полярных координатах:
Будем решать задачу Дирихле для уравнения Пуассона.
в прямоугольной области , .
Для численного решения данной задачи применим метод SOR (метод последовательной верхней релаксации). Вначале используем метод конечных разностей. Для этого разобьем отрезок [a,b] на K равных интервалов длиной , а отрезок [c,d] - на L интервалов той же длины. Пусть при этом , . Тогда
Таким образом, мы разбили весь прямоугольник , на квадраты со стороной ,получив при этом сетку (решетку). Будем искать нужную нам функцию в узлах этой решетки методом конечных разностей. То есть будем искать эту функцию в виде матрицы
Запишем вначале в одномерном случае рамках метода конечных разностей производную функции U по х в точке . Получим
Для второй производной в точке получим
Аналогичным образом, для второй производной в точке имеем
Перейдем теперь к двумерной области и найдем вторые производные в точке , .
,
Подставив полученные таким образом вторые производные в уравнение Пуассона, имеем
Коэффициенты перед матричными элементами в данном случае равны 1 и 4. Однако в общем случае (полярные координаты, например) мы должны записать
и вычислить A, B, C, D, E.
Перепишем полученное нами уравнение в следующем виде:
Разность между левой и правой частями уравнения называется невязкой
Или
Подставляя в правую часть (3) уравнение (1) имеем
Уравнение (4) является тождеством. Его можно использовать в случае, когда решение ищется методом последовательных приближений
Это и есть основное уравнение метода SOR. Параметр обычно выбирается не тождественно равным 1, изменяющимся параметром от 1 до 1.5. Это делается для лучшей сходимости метода. В качестве начального приближения берется функция, тождественно равная 0 внутри области определения и равная граничным условиям на границе области.
Будем искать решение этой задачи методом Кранка-Николсона. Рассмотрим вначале одномерное нестационарное уравнение диффузии.
Запишем производную по времени в k-й точке х в n-й момент времени в виде:
Здесь - шаг по времени. В правой части вторую производную по координате можно взять как в момент времени ,
так и в момент времени
Суть метода Кранка-Николсона заключается в том, что производная в правой части берется как среднее арифметическое от производных в точках и . То есть
Для двумерного уравнения диффузии
в рамках метода Кранка-Николсона запишем вторые производные по х и по у
как среднее арифметическое от производных в точках и ,
перенесем все значения функции в момент времени влево, а в -й вправо. Получим.
,
Или
(1)
(2)
Схема расчета по методу Кранка-Николсона такова:
МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ
Базисные функции метода конечных элементов
Суть метода в том, что искомая функция ищется интерполяцией с помощью базисных функций. В этом случае отпадает необходимость расчета производных методом конечных разностей.
Пусть отрезок [a,b], на котором определяется искомая функция, разбит N точками на (N-1) равных отрезков (конечных элементов) длиной . На рисунке показан случай, когда N = 4.
Значения функции в этих точках равны , , .... , . на каждом из конечных элементов функцию u(x) можно представить с помощью базисных функций :
|
1 элемент |
|
|
2 элемент |
|
|
(N-1) элемент |
Суммируя по всем элементам, получим:
На каждом отрезке изменяется от 0 до 1.
Вид функций для N = 4 приведен ниже.
(1)
Рассмотрим перечисленные выше шаги.
- весовая функция.
(2)
Пусть g есть функция только х. Тогда в случае, когда область интегрирования - отрезок [a,b], интегральное уравнение имеет вид:
(3)
(4)
Тогда
(5)
В этом случае интегральное уравнение (3) преобразуется к виду
(6)
Или
(7)
Интегралы по элементам
Представим u(x) в виде
Используем аппроксимацию Галеркина:
Теперь уравнение (7) записывается как:
(8)
Здесь u - вектор значений температуры в узлах, элементы матрицы Е (Е - матрица жесткости) можно вычислить по формуле:
Разделим область [a,b] на (N-1) элементов равной длины. Пусть a =1, b=1, N = 3.
Тогда для каждого элемента
Тогда
Рассчитаем E11, E12, E21, E22. Теперь необходимо путем ансамблирования рассчитать глобальные значения элементов матрицы Е. Получим:
Рассмотрим теперь правую часть уравнения (8)
(9)
В случае, если заданы потоки тепла на концах стержня и , элементы вектора достаточно просто.
Рассмотрим случай , когда в качестве граничных условий заданы значения температуры. Тогда
Аналогичным образом вычисляется . Теперь необходимо перенести известные значения un из левой части системы (8) в правую, а неизвестные - из правой части в левую. При этом система уравнений (8) будет переопределена и ее необходимо решать минимизируя функционал Q.
Двумерные базисные функции
Выбрав таким образом двумерные базисные функции, искомую функцию на квадратном конечном элементе можно представить как:
Для реализации метода конечных элементов нам необходимо:
1 шаг.
2 шаг.