Методы решения дифференциальных уравнений гидродинамики Текст научной статьи по специальности «Математика»

УДК 532.5:517.9 ББК 22.253.3 М-54

Меретуков Заур Айдамирович, кандидат технических наук, докторант кафедры технологии, машин и оборудования пищевых производств Майкопского государственного технологического университета, e-mail: zamer@radnet.ru;

Заславец Александр Алексеевич, соискатель кафедры машин и аппаратов пищевых производств Кубанского государственного технологического университета, т.: (861)2752279;

Кошевой Евгений Пантелеевич, доктор технических наук, профессор, заведующий кафедрой машины и аппараты пищевых производств Кубанского государственного технологического университета, e-mail: Koshevoi@kubstu.ru;

Косачев Вячеслав Степанович, доктор технических наук, профессор кафедры машин и аппаратов пищевых производств факультета машиностроения и автосервиса Кубанского государственного технологического университет, т.: (861) 275-22-79.

МЕТОДЫ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ГИДРОДИНАМИКИ

(рецензирована)

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

Ключевые слова: непрерывная функция, функция Хэвисайда, функция округления, метод Галеркина.

Meretukov Zaur Aydamirovich, Candidate of Technical Sciences, doctoral student of the Department of Technology, Machinery and Equipment for Food Industries Maikop State Technological University, e-mail: zamer@radnet.ru;

Zaslavets Alexander Alexeevich, seeker of the Department of Machines and Apparatus for Food Production, Kuban State Technological University, tel: (861) 2752279.

Koshevoi Eugene Panteleevich, Doctor of Technical Sciences, professor, Head of the Department of Machines and Apparatus of Food Production, Kuban State Technological University, e-mail: Koshevoi@kubstu.ru;

Kosachev Vyacheslav Stepanovich, Doctor of Technical Sciences, professor, professor of the Department of Machines and Apparatus of Food Production, Kuban State Technological University, tel.: (861) 2752279.

METHODS OF SOLVING DIFFERENTIAL EQUATIONS OF HYDRODYNAMICS

(reviewed)

The aim is to obtain the simplified method for solving differential equations of hydrodynamics with a high accuracy of the results that the use of a regression function will determine the required number offinite element of the trial function.

Key words: continuous function, the Heaviside function, the function of rounding, the Galerkin method.

Одной из важных задач в анализе континуальных систем, например в задачах гидродинамики [1], является преобразование непрерывных функций в набор конечных элементов [2]. Метод конечных элементов позволяет строить удобную схему алгебраических уравнений относительно узловых значений искомой функции. При этом приближенная аппроксимация решения осуществляется на типовом полиномиальном элементе. В данной работе предлагается дополнить используемые полиномиальные элементы ступенчатой функцией Хэвисайда и функциями округления, что позволяет не только упростить и формализовать запись пробной кусочно-непрерывной функции, но и избежать интегрирования при ее использовании для решения континуальных задач методом Галеркина [3]. Специфическим в методе конечных элементов является построение семейства функций, определяемых ограниченным числом параметров. Например, семейство функций u(x) при Xmin < х < Xmax. Интервал Xmin ... Xmax представляет собой одномерную область существования решения задачи, которая

разбивается на конечное число частей (элементов), соединяющихся между собой и с концами интервала в узловых точках (узлах) XI. В пределах каждого элемента задается функция виде линейного полинома. Она определяется своими значениями u(Xi) в узлах и на концах элемента. Учитывая, что в континуальной задаче функция является непрерывной, то ее значения в каждом узле для соседних элементов должны совпадать. Для этого введем функции округления: _х} - функция пола, которая определяется как наибольшее целое, меньшее или равное x, а именно _xJ = п о x - 1 < п < x; Гх1 -функция потолка, которая определяется как наименьшее целое, большее или равное x, а именно _xJ = п о x < п < x +1. Используя эти функции округления, введем семейство кусочно-линейных непрерывных функций следующего вида:

г х — а

и(х) = {_Ф(х — а )} — ГФ(* — Ь)1}

с + (<^ — с) ■

Ь — а

(1)

где Ф(х) - функция Хэвисайда, единичная ступенчатая функция, чьё значение равно нулю для отрицательных аргументов и единице для положительных аргументов и 1/2 в нулевой точке; a, Ь -интервал окна функции u(x) на котором она отлична от нуля (а < Ь); с, й - параметры уравнения прямой на интервале а, Ь. Функция и(х) на границах интервала (а, Ь) обращается в ноль. Используем (1) для описания пробной функции на произвольном интервале, задаваемом полушириной И = (Ь -а)/2 одномерного единичного конечного элемента относительно узла X В этом случае параметры кусочно-линейной непрерывной функции (1) соответственно равны а = X] - И, Ь = X] + И, с = 0, й = 1, а функция приобретает вид:

<Рг (Х) = 1Ф[Х — (Хг — И)]} — ГФ(Х — Х1 )1}- Х + И Х' +

+ ■

{ГФ(х — X{ )1 — _Ф[х — (Хг + И)]}}

И

И — х + X ■

(2)

И

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

Метод конечных элементов оперирует в качестве ф; (х) кусочно-полиномиальными функциями, отличными от нуля в пределах небольшого числа элементов вблизи узла X Поскольку и(х) по своему физическому смыслу должна быть непрерывной функцией, выберем ф(х) в виде кусочно-линейных функций, отличных от нуля на двух соседних элементах относительно выбранного (Рис. 1). Каждая такая функция ф/х), г = 0, 1, 2, ..., п-1, равна единице в X] и нулю во всех остальных узлах. При этом набор функций и(х) будет состоять из непрерывных функций, линейных в пределах элементов с изломами в узлах и определяемых своими узловыми значениями иг-. Каждую из таких функций можно изобразить в виде ломаной линии. Метод конечных элементов заменяет задачу отыскания функции на задачу отыскания конечного числа ее приближенных значений в отдельных точках-узлах.

X

Рис. 1. Семейство кусочно-линейных функций (2) для Х0 = (-1/2), Хг = 0, Х2 = (-1/2) с полушириной конечного элемента Л = 1/2 используемых при решении континуальной одномерной задачи методом конечных элементов на отрезке [-1; 1]

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

С уменьшением максимального размера элементов увеличивается число узлов и неизвестных узловых параметров. Вместе с этим повышается возможность более точно удовлетворить уравнениям задачи и тем самым приблизиться к искомой функции. Рассмотрим применение метода для описания одномерного профиля. В этом случае пробная функция определяемая уравнением (3) представлена линейной комбинацией функций (2) с коэффициентами щ:

п—1

и(х) = 2 и -ф М (3)

г=0

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

п—1

X и 'Фг (Х)= 1 (4)

1=0

и умножая левую и правую часть равенства (4) на ф(х) гдеу=0, 1, 2, ..., п-1:

п

р. (х )•£ ир (х ) = 1-р. (х)

(5)

і=1

интегрируем систему уравнений (5) по области существования решения Хтіп ... Хтах, согласно методу Галеркина получаем систему алгебраических уравнений, используемую для определения коэффициентов П{.

х тах п—1 тах

| Р і (х) - ^ и ' Р (х) іїх =| р (х ')dx

Xтіп I- >=0 J Хт,п

(6)

При этом определенные интегралы имеют простые алгебраические выражения во всех узлах X, на интервале существования решения:

°, / < j — 1

к • • 1

—, г = j — 1

6 2

dx = ^ — • к, г = ] (7)

3

к • • ,

—, г = ] +1 6

о, г > j +1

х„

X т

г П—1

I Р. (х )'Х иі Р (х)

і=0

х„

| Р (х )dx = к

(8)

х„

где к = (Хтіп - Хтах)/(п + 1) - полуширина конечного элемента. Используя соотношения (7) и (8) в системе уравнений (6) получаем слева трех диагональную разреженную матрицу коэффициентов умноженную на вектор неизвестны {и}, а справа вектор из элементов к.

(9)

1 к к 6 0 0 0 и0 к

к 6 1 к к 6 0 0 и к

0 к 6 1 к к 6 0 - и2 = к

0 0 к 6 І к к 6 ип—2 к

0 0 0 к 6 1 к ип—1 к

Умножая левую и правую часть матричного уравнения (9) на (1/к) получаем матричное уравнение не использующее параметр (к).

2 3 1 6 0 0 0 Щ 1

1 6 2 3 1 6 0 0 и 2 1

0 1 6 2 3 1 6 0 - и3 = 1

0 0 1 6 2 3 1 6 ип—1 1

0 0 0 1 6 2 3 ип 1

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

(11)

2 3 1 6 0 0 0 Щ и (X 0)

1 6 2 3 1 6 0 0 и 2 и (X1,)

0 1 6 2 3 1 6 0 - и3 = и (X 2 )

0 0 1 6 2 3 1 6 ип—1 и (Xn—2 )

0 0 0 1 6 2 3 ип и (Xn—,)

Здесь матрица (11) является трех диагональной, и решение может быть получено методом Гаусса, матричным методом (умножением на обратную матрицу) или прогонки. Для систем малой размерности метод решения не существенен, но с увеличением точности решения необходимо увеличивать число конечных элементов на области существования решения, а это приводит к увеличению числа неизвестных в уравнении (11). Если требуется большая точность решения необходимо использовать метод прогонки, позволяющий решать системы большой размерности без существенных ошибок округления. Рассмотрим применение этого метода в данном случае. Метод прогонки является двух шаговым. Вначале вычисляем вспомогательные величины а,;, зависящие от элементов трех диагональной матрицы из правой части (11):

11 1

Сп —--------------, ссл —

0 4 1

4 + сс

•,Сп—1 =

4 + с

(12)

"0 ' 1 ^п—2

Рекуррентные соотношения (12) фактически представляют собой цепную дробь следующего

вида:

1

4 —

(13)

4 — — 4

Следовательно (13) может быть представлена рекуррентной формулой образованной следующей числовой последовательностью:

{Яг—1 )2 — 1

Чг = —■

Чг —2

(14)

где до = 4; ^ = 15.

В этом случае элементы вектора (12) задаются следующими соотношениями:

Яо,...,а=— ^

Я ‘ Яг

с0 = —

1

Ч0

<С — —

(15)

Следующим этапом первого шага метода прогонки является расчет коэффициентов ";, которые для уравнения (1о) (единичный прямоугольный профиль) определяются следующими соотношениями:

3 Д = 6 — ро д = 6 — Д—1

2 , Д ................ "

4 + с

4 + с

(16)

-0 ^ ' ^/—1 Используя соотношения Эйлера преобразуем (16) к удобному для вычислений рекуррентному виду образованному следующими числовыми рядами:

=

1

а0 = 3, ах = 6, ..., а = Ьо = 2, Ь1 = 5 •••, Ь =

з -(-1)' 2

3-(-1)'

2

а'-1 + а'-2

• Ь'-1 + Ь'-

2

(17)

(18)

Используя (17) и (18) получаем рекуррентную формулу для расчета элементов вектора (16):

а,

— ___о н — _______1 к —

,0 V"”" Ь (19)

Определив из (15) и (19) значения вспомогательных векторов переходим ко второму шагу метода прогонки - расчету весовых коэффициентов иг- по формулам:

6 — Дп—2 г>

ип-1 =--, ип—2 = ™п—2 • ип—1 + Дп—2 =

4 + ^п—2

ио =ао • и1 + До (20)

Используя значения весовых коэффициентов иь для пробной функции (3) получаем аппроксимацию этой функцией исходного единичного прямоугольного профиля:

п—1

(X ) = /\и, • (рк (х, п)]

п — 1

¥ (х) = / \ы

/Л11 к

к=0

(21)

где ф(х, п) определяется следующей формулой:

Рк (хп)

2•к 1

+ Ф х---------+ 1

V п +1 у

Интегрирование функции числовую последовательность:

Г 2 • к + 2 Л Г 2 • к + 4 Л { -> Л , п х п • х 3

Ф х +1 — Ф х +1 к +

V п+1 ) V п+1 ) (М 2 2 2

+

— Ф

х — ■

2 • к + 2 ,

-------+1

п + 1 у

п X п • X 1

— к + — +--+ —

V 2 2 2 2у

(22)

(21) на интервале прямоугольного профиля дает следующую

„ 4 г 6 „ 84 47 4 56 9 0 248 3411 103 0 8 1286

£ :=- { ;=_ г :=- ^ :=- { := — ^ :=- ^ :=_ ^ :=— { •=----------------------------------------------- { :=

2 5 3 7 4 95 5 52 6 497 7 97 8 265 9 3620 10 10879 11 1351

— (23)

Таким образом, предложенный конечный элемент (22) для аппроксимации континуальных функций позволяет использовать числовые ряды для оценки точности описания континуальной функции, не прибегая к численному интегрированию. Рассмотрим данную числовую последовательность (23) в виде графика в осях (1-:Р(п) -1/п} - (рисунок 2).

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

18%

16%

14%

12%

10%

8%

6%

4%

2%

0%

1,07691 x *-0,4168 R2 = 2xs - 0,39 0,99998 B94x2 + 0

У Ou wU 1 X

0 0,06 0,1 0,15 0,2 0,26 0.3 0.36 0,4 0,45 0.5

1 /п

Рис. 2. Экстраполяция числового ряда невязок - [1-/(п)] от обратного числа

конечных элементов - (1/п)

Литература:

1. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: пер. с англ. В 2 т. Т.1. М.: Мир, 1990. 384 с.

2. Зенкевич О., Морган К. Конечные элементы и аппроксимация: пер. с англ. М., 1986. 318 с.

3. Флетчер К. Численные методы на основе метода Галеркина: пер. с англ. М.: Мир, 1988. 352 с.

References:

1. Anderson D., Tannehill J., Pletcher R. Calculating hydroechanics and heat transfer: trans. from English. In 2 vol. Vol.1. M.: Mir, 1990. 384 p.

2.Zenkevich O., Morgan K. Finite elements and approximation: trans. from English. M., 1986. 318 p.

3. Fletcher K. Numerical methods based on the Galerkin method: trans. from English. M.: Mir, 1988.

352 p.