CC BY
18
Bemmx&TWHT/Proceedmgs of VSUET, Т. 80, № 2, 2018-
Оригинальная статья/Original article_
УДК 519.635.1
DOI: http://doi.org/10.20914/2310-1202-2018-2-114-118
Интегрирование бигармонического уравнения по неявной схеме
Михаил И. Попов 1 mihail semilov@mail.ru
1 Воронежский государственный университет инженерных технологий, пр-т Революции, 19, г. Воронеж, 394036, Россия Реферат. В статье представлено пошаговое построение конечно-разностной схемы для неоднородного бигармонического уравнения при нулевых граничных условиях, наложенных на искомую функцию и ее частные производные первого порядка. Конечно-разностная схема основана на квадратном двадцатипятиточечном шаблоне и имеет неявный характер. На равномерной сетке с помощью разложения функции в ряд Тейлора с остаточным членом в форме Лагранжа вычислена погрешность аппроксимации бигармонического оператора разностным аналогом и погрешность аппроксимации граничных условий, наложенных на частные производные первого порядка. Граничные условия, наложенные на искомую функцию, выполняются точно. Конечно-разностная схема аппроксимирует краевую задачу со вторым порядком точности по шагу сетки. С помощью системы компьютерной алгебры Maple получено решения задачи для различных шагов сетки. Выявлена зависимость минимума функции и времени расчета от числа значимых цифр. Найдено оптимальное число значащих цифр. Проведен анализ скорости сходимости численной схемы. Установлена зависимость минимального значения функции и времени расчета от величины шага сетки. Найдено оптимальное значение шага. Построены трехмерный график решения и его профили в серединных сечениях.Указаны преимущества разработанной конечно-разностной схемы. Полученные результаты отвечают физическому смыслу задачи и согласуются аналогичными численными и приближенно-аналитическими решениями.
Ключевые слова
^раевая^адача^бигармоническое^равнение^
The integration of a biharmonic equation by an implicit scheme
Mikhail I. Popov 1 mihail_semilov@mail.ru
1 Voronezh state university of engineering technologies, Revolution Av., 19 Voronezh, 394036, Russia__
Summary.The paper presents a step-by-step construction of a finite-difference scheme for a heterogeneous biharmonic equation under zero boundary conditions superimposed on the desired function and its first-order partial derivatives. The finite-difference scheme is based on a
square twenty-five-point pattern and has an implicit character. On analytical grid, the error of approximation of the biharmonic operator by the difference analog and the error of approximation of boundary conditions imposed on the first order partial derivatives are calculated by the expansion of the function in the Taylor series with the remainder term in the form of a Lagrange. The boundary conditions imposed on the
sought function are satisfied precisely. A finite-difference scheme approximates a boundary value problem with a second order of accuracy over the mesh step. With the help of the Maple computer algebra system the solutions of the problem for different grid steps are obtained. The dependence of the minimum function and calculation time on the number of significant digits is revealed. The optimal number of significant digits is found. The convergence rate of the numerical scheme is analyzed. The dependence of the minimum value of the function and the calculation time on the value of the grid step is established. The optimal step value is found. A three-dimensional graph of the solution and its profiles in the middle sections are constructed. The advantages of the developed finite-difference scheme are indicated. Obtained results correspond to the physical meaning of the problem and are consistent with similar numerical and approximate analytical solutions. Keywords: boundary value problem, biharmonic equation, finite-difference scheme, approximation error
Введение
Краевые задачи для бигармонического уравнения возникают в теории упругости при моделировании изгибов тонких пластин [1-2], а также в гидродинамике при моделировании свободно-конвективных течений [3-5].
Наиболее известное решение неоднородного бигармонического уравнения, так называемого уравнения Софи Жермен, получены Анри Навье в виде двойного тригонометрического ряда. Однако оно справедливо лишь для пластинок, шарнирно опертых по контуру. Ряды в этом решении сходятся медленно. Более общее решение предложил Морис Леви, но и им не исчерпываются все постановки граничных условий.
В связи с отсутствием аналитических решений широкое развитие получили численные методы интегрирования краевых задач для бигар-монического уравнения. Наряду с классическими [6-7] используются и современные методы, такие как метод конечного множества точек [8-9].
Для цитирования Попов М.И. Интегрирование бигармонического уравнения по неявной схеме // Вестник ВГУИТ. 2018. Т. 80. № 2. С. 114-118. ао1: 10.20914/2310-1202-2018-2-114-118
Построение конечно-разностной схемы
Рассмотрим краевую задачу для неоднородного бигармонического уравнения с граничными условиями первого и второго рода в квадратной области В = [0;1] х [0;1].
д4Ф(х,у) | 2д4Ф(х,у) | д4Ф(х,у) = ^ дХ4 дх2д¥2 дУ4 ' Ф(0,7) = Ф(1,У) = Ф(х,0) = Ф(Х,1) = 0, (2)
дФ(0,Г) _ дФ(1,Г) _ дФ(X,0) _ 5Ф(X,1) _ Q (3)
5X 5X 5Y 5Y
Для построения конечно-разностной схемы заменим область D решения системы сеткой с шагами дх , a y ,
Wx,AY = {(X,Yj) = (i AX, jAY)}, / = 0, .., n j =0, ,.,m
Вместо функции непрерывного аргумента на D рассматривается функция дискретного аргумента Ф(X,Y дY, которая обозначается
For citation
Popov M.I. The integration of a biharmonic equation by an implicit scheme. Vestnik VGUIT [Proceedings of VSUET]. 2018. vol. 80. no. 2. pp. 114-118. (in Russian). doi:10.20914/2310-1202-2018-2-114-118
Ветшк&ТУИТ/Ргоахб^ о/ Т. 80, № 2, 2018-
Выбор разбиения осуществлен таким образом, чтобы в него попадали точки границы области и центр области.
Для построения конечно-разностного аналога частных производных д4Ф(Х,У)/дХ4 и Э4Ф(Х, У)/дУ4 используем центрально-разностный оператор 2-го порядка примененный дважды по соответствующей переменной:
Конечно-разностный аналог смешанной производной
д 4Ф
дХ 2дУ2
144(ДХ )2(АУ )2
[(Ф,+ 2,1+2 -16Ф,+2,1+1 +
30ф;+27 -16фг+2,^_1 + Ф,+ад_2
'¡+1,1+2 -16Ф,+1,1+1 +
д4ф дХ4
. Ф,+2,1 - 4Ф,+1,1 + 6Ф,, 1 - 4Ф,-1,1 + Ф,-2,1 (ДХ)4
(4)
)-16 к
+30ФМ,1 - 16Ф;+1,1 -1 +Фм,1-2 )+ 30 (Ф,1+2 - 16Ф.1 + 1 + (6) +30Ф;,1 - 16Фц +Ф;,1-2 )- 16 (Фм,1+2 -16Ф;-1,1-2 + ) + (Ф,-2,1+2 - 16Ф,-2,1+1 + -16Ф* , + Ф*
+30Ф ,-1,1 - 16Ф,-1,1 -1 +ФМ,1 -2
+30Ф*
д 4Ф дУ4
_ Ф/, 1+2 - 4Ф,, 1+1 + 6Ф, 1 - 4Ф,, 1-! +Ф, 1-2 (5) 1 (АУ)4
'¡-2,1 16Ф,-2,1-1 1 Ф,-2,1 -2 Подставив в уравнение (1) аппроксимации производных (4)-(6), получим разностное уравнение
Ф,+2,1 - 4Ф,+1,1 + 6Ф,,1 - 4Ф,-1,1 + Ф,-2,1 | Ф,, 1+2 - 4Ф,,1+1 + 6Ф,,1 - 4Ф,, 1 -1 + Ф,, 1 -2 +
(ДХ)4
"[(Ф,+2,1+2 - 16Ф,+2,1+1 + 30Ф,+ 2,1 16Ф,+2,1-1 + Ф,+2,1-2
(АУ)
)-16 (Ф,+1,1+2 -16Ф
¡+1,1+1 +30Ф,+1,1 -16Ф,+1,1-1 + Ф,+1,1 -2 ) +
72(ДХ )2(Д 7)
+30 (Ф,,1+2 - 16Ф,,1+1 +30Ф,,1 -16Ф,, 1-1 +Ф,,1-2 )- 16 (Ф,-1,1+2 16Ф,1,1 -2 +30Ф^1,1 - 16Ф^,1-1 +Ф^1,1-2 ) +
+ (Ф^-2,1+2 -16Ф?-2,1+1 +30Ф^-2,1 16Ф^-2,1 -1 +ФЦ1 -2 )] = -1
Граничные условия (2) накладываемые Тогда уравнение (7) примет вид
на искомую функцию приобретут вид
Ф = Ф = Ф = Ф =0, , = 0,п, 1 = 0,п. (8)
0,1 п,1 1,0 ¡,п ' ' ' ^ ' V"/
Граничные условия накладываемые на частные производные первого порядка получим, используя разностные операторы для левой и правой границ
_1_ /1гТ» _ гТ»
, И -1 (9)
) + (7)
дФ -3 Ф, 0 + 4 Ф,д -Ф,,2
дХ 0 2АХ
дФ 3Ф. ¡,п -4Ф. , , ,п-1 + Ф,,п-2
дХ ¡,п 2АХ
I = 1, , = 1.....п-1 (10)
Для верхней и нижней границ будем иметь соответственно
дФ д У дФ 57
_ -3Ф0,1 + 4Фю-Ф2,]
0.1
2ДХ
3Ф - 4Ф , +Ф „
_ 3Фп, 1 4Фп-1,1 +Фп-2,1
1 2ДХ
1 = 1,. 1 = 1,
, п -1 (11) .., п -1(12)
Я*Ф = -1
Чтобы определить погрешность аппроксимации образуем разность г = и - V , где и -решение задачи (7), (8), (13), (14), а V - решение задачи (1)-(3). Подставляя и = г + V в уравнение, получим для г задачу
Ви = В(г + V) = Вг + Bv = -1, Вг = -ев
во внутренних узлах сетки и для г выполнены граничные условия (8), (13), (14), где ев = В,у + 1 -погрешность аппроксимации задачи (1)-(3) схемой (7), (8), (13), (14). Так как Д\>+1 = 0, то
ев = Bv +1 - А2 — 1 = В^- А2
Вычислим погрешность аппроксимации почленно, используя разложение по формуле Тейлора с остаточным членом в форме Ла-гранжа, в окрестности ,, 1 -того узла:
Заменим в уравнении (3) производные их конечно-разностными аналогами (9) - (12). Тогда, учитывая уравнение (8), граничные условия на производные примут вид
, = 1,...,и -1 (13) 1, 1 = 1,...,п -1 (14)
д4 V
Ф,
-Ф,
Ф , = —Ф
,п 1 4 •
Ф„
. =-Ф„
дХ4
к! _д!г
6 дХ
д4 V
д у 4
д4 v( Х ,У)
1
дХ4
6[к(, + 7), ]к], 0 < 77 <1,
•, 1
д4 v( Х ,У) дУ4
В итоге непрерывная краевая задача (1) - (3) свелась к конечно-разностной схеме (7), (8), (13), (14).
Вычисление погрешности аппроксимации
На равномерной сетке ДХ = ДУ = к обозначим
к!
6 ду6
[¡к, к( 1 + а)], 0 < а <1,
д4 v( Х, У)
дХ2 дУ2
дХ 2дУ2
В =
д4Ф дХ4
+ 2-
д 4Ф дХ 2дУ2
д4Ф "дУ4
1
1
+
17 h J 58v r , 58vr.,;/. ч1,
-i—¿tUi + n)h, hj 1+—-\ih, h( j +
756 I 578 1 J 578 L J|
58 v
34h4 _
189 W757
58v
[(/ + n)h, h(j + ст)] +
+
5775X 167h4 I 58v
270 l5X6572
58v
[(/ + n)h, h(j + CT)]|-
[(i + rî)h, h(j + ст)] +
+
576 5X2 34h4 I 58 v
27 15X5 573
58v
[(i + rj)h, h(j + CT)]|-[(i + n)h, h( j +ст), кт~] +
+
575 5X3
[(i + n)h,h(j + CT),кг] r- +
85h 58v
54 5X4574 Обозначим
5 6v( X ,7 )
[(i + n)h,h(j + ct)],0 <n < 1,0 <a < 1.
M, = max
1 (X ,7 )eD
M = min_
3 (X ,7 )eD
M. = max_
4 (X ,7 )eD
5X6 5m+nv( X ,7 )
д 6v( X ,7 )
5Xm 57 n 5m+nv( X ,7 )
5Xm 57n
Mn = max ,
(X,7)eD 576
, m, n = 1,3,5,7, m + n = 8; m, n = 2,4,6,8, m + n = 8.
Поскольку бигармонический оператор симметричен М1=М^. Обозначим Ml /w /з . Л/2 = (5398/945Ж4 -(1088/189)М"3, тогда
I SB \<Mjf +MnhA.
(15)
Таким образом, разностный оператор аппроксимирует бигармонический со 2 порядком.
Граничное условие (2) на искомую функцию с учетом (8) аппроксимируется точно. Погрешность аппроксимации граничного условия (3) уравнениями (13) и (14) оценим используя разложение по формуле Тейлора с остаточным членом в форме Лагранжа, в окрестности /', у -того узла: ду
дХ у"
Обозначим
д3у( X ,У )
-— = -h -^-V[h(i + n), jh], 0< n <1.
5X 3 5X3
M = max_
( X ,7 )eD
5X3
тогда I ет I< Мк2.
Итак, граничное условие (2) аппроксимируется также со вторым порядком, поэтому конечно-разностная схема (7), (8), (13), (14) аппроксимирует краевую задачу (1)-(3) со вторым порядком точности.
Анализ решений
Решение задачи (7), (8), (13), (14) вычислено с помощью системы компьютерной алгебры Maple для равномерной сетки с шагом h = 1120. На рисунке 1 показан график решения, а на рисунке 2 его профиль в сечениях плоскостями X = 0.5 и Y = 0.5 . На рисунке 3 показаны профили решения в сечениях плоскостью Y = X и Y =1 - X . По графикам видно, что решение симметрично относительно центральной точки (0.5;0.5), в которой достигается минимум равный -1.22 • 10-3.
В таблице 1 показана зависимость минимума функции от числа значимых цифр в расчете, т.е. фактически от погрешности округления. Расчет выполнен для сетки размерностью 81 х 81 точек. В столбце относительная разница минимумов функции показано насколько процентов результат отличается от предыдущего.
Здесь же приведено время расчета и егоот-носительное изменение. По таблице видно, что при увеличении числа значащих цифр с 14 до 16 время расчета увеличивается на 166% а точность всего на 3 • 10-6 %. Эти данные говорят о том, что дальнейшее уменьшение погрешности округления приведет к значительному увеличению времени расчета при неизменном результате. Таким образом, 14 значащих цифр - оптимальным значение.
В таблице 2 отражена зависимость минимального значения функции от величины шага сетки и относительная разница минимумов функции. Здесь же приведено время расчета и его относительное изменение.Расчет выполнен с 14 значащими цифрами. Уменьшение относительной разницы минимумов функции по экспоненциальному закону свидетельствует о быстрой сходимости конечно-разностной схемы.
Рисунок 1. График решения Figure 1. Solution graph
Рисунок 2. Профиль решения в сечении X = 0.5 Figure 2. Solution profile in section X = 0.5
Уменьшение относительного приращения времени расчета позволяет сделать расчеты более точными. Однако уменьшение шага сетки потребует увеличения значащих цифр, поэтому значение шага 1/120 является оптимальным
Рисунок 3. Профили решения в сечениях y = х (красный) и y = 1 — х (синий) Figure 3. Solution profile in section y = x(red) and y = i-х (blue)
Полученный результат согласуется с численным решением данной задачи по явной схеме [10] и приближенно-аналитическим решением [11]. Значение функции в центральной точке для которых равны —1.05-10_3 и -1.25 -10 3 соответственно.
Погрешность округления Rounding error
Таблица 1. Table 1.
Число значащих цифр Number of significant digits Минимум функции The minimum of a function Относительная разница минимумов функции, % The minimum of the relative difference in the minimum of a function, % Время расчета Calculation time Относительная разница времени расчета, % Relative time difference calculation, %
8 -0.0012409828 333.687
10 -0.001202252713 -3.12 324.56 -2.74
12 -0.00120187713348 -0.03 326.96 0.74
14 -0.0012018733801815 -0.0003 332.969 1.84
16 -0.001201873342732467 -3-10-6 888.996 166.99
Таблица 2.
Скорость сходимости
Table 2.
Rate of convergence
Шаг сетки Grid spacing Минимум функции The minimum of a function Относительная разница минимумов функции, % The minimum of the relative difference in the minimum of a function, % Время расчета Calculation time Относительная разница времени расчета, % Relative time difference calculation, %
1/20 -0.0010105920924147 1.357
1/40 -0.0011381686297206 12.62393987 14.165 943.8467208
1/60 -0.0011806597334856 3.733287131 86.44 510.2364984
1/80 -0.0012018733801815 1.796762133 332.969 285.2024527
1/100 -0.0012145881472324 1.057912357 1090.431 227.4872436
1/120 -0.0012230583683455 0.6973739223 3259.407 198.9099723
Заключение
С помощью неявной конечно-разностной схемы получено решение неоднородного бигармонического уравнения. Путем разложения функции в ряд Тейлора вычислена погрешность аппроксимации бигармонического оператора разностным аналогом и погрешность аппроксимации граничных условий. Проведен анализ скорости сходимости схемы. Преимущества использования данной неявной схемы по сравнению с явной заключаются в ее абсолютной ЛИТЕРАТУРА
1 Завьялов В.Н., Мартынов Е.А., Романовский В.М. Основы строительной механики пластин: учебное пособие. Омск: СибАДИ, 2012. 116 с.
2 Шафарец Е.Б., Шафарец Б.П. Свободная конвекция учет некоторых физических особенностей при моделировании конвективных течений с помощью вычислительных пакетов // Научное приборостроение. 2014. Т. 24. №2. С. 43-51.
3 Гоц А.Н. Численные методы расчета в энергомашиностроении. Владимир: Изд-во ВлГУ, 2013.182 с.
4 Jani S., Mahmoodi М., Amini М., Jam J. Numerical investigation of natural convection heat transfer in a symmetrically cooled square cavity with a thin fin on its bottom wall // Thermal science. 2014. V. 18. №. 4. P. 1119-1132.
5 Gros Т., Revnic C., Pop I., Ingham D.B. Free convection heat transfer in a square cavity filled with a porous medium saturated by a nanofluid // International Journal of Heat and Mass Transfer. 2015. V. 87. P. 36^11.
бАлгазин С.Д. Численные алгоритмы классической математической физики. М.: Диалог-МИФИ, 2010. 240 с.
7 Mu L., Wang J., Ye X. Effective implementation of the weak Galerkin finite element methods for the biharmonic equation // Computers & Mathematics with Applications. 2017. V. 74. №. 6.P. 1215-1222.
8 Doss L. J. Т., Kousalya N. Finite Pointset Method for biharmonic equations // Computers & Mathematics with Applications. 2018. V. 75. №. 10. P. 3756-3785.
9 Doss L. J. Т., Kousalya N, Sundar S. A Finite Pointset Method for Biharmonic Equation Based on Mixed Formulation // International Journal of Computational Methods. 2017. P. 1850068.
10 РяжскихВ.И, Слюсарев М.И., Попов М.И. Численное интегрирование бигармонического уравнения в квадратной области // Вестник Санкт-Петербургского университета. 2013. № 10. V. 1. Р. 52-62.
11 Попов М.И., Соболева Е.А. Приближенное аналитическое решение внутренней задачи кондуктивно-ламинарной свободной конвекции // Вестник ВГУИТ. 2016. № 4. С. 78-84.
устойчивости. Доказательство устойчивости явных схем порой весьма затруднительно. Отпадает необходимость вычислять нормы разностных операторов для определения оптимального для скорости сходимости шага. Объем вычислений существенно меньше, так как значения в узлах вычисляются один раз, а не на каждом шаге. Это позволяет экономить ресурсы процессора и время расчета, а также уменьшить погрешность результата за счет уменьшения числа округлений.
REFERENCES
1 Zavialov V.N., Martinov E.A., Romanovskyi V.M. Os-novi stroitelnoy mehaniki plastin [Fundamentals of structural mechanics ofplates]. Omsk, SibADI, 2012, 116 p. (in Russian)
2 Shafarets E.B., Shafarets B.P. Free convection taking into account some physical features when modeling con-vective flows using computational packages. Nauchnoe pribo-rostroenie [Scientific instrument engineering], 2014, vol. 24, no 2, pp. 43-51.(in Russian)
3 Gots A.N. Chislennie metodi raschota v ener-gomashinostroenii [Numerical methods of calculation in power engineering]. Vladimir, VlGU, 2013,182 p.
4 Jani S., Mahmoodi M., Amini M., Jam J. Numerical investigation of natural convection heat transfer in a symmetrically cooled square cavity with a thin fin on its bottom wall. Thermal science, 2014, vol. 18, no. 4, pp. 1119-1132
5 Gros T., Revnic C., Pop I., Ingham D.B. Free convection heat transfer in a square cavity filled with a porous medium saturated by a nanofluid. International Journal of Heat and Mass Transfer, 2015. vol. 87. pp. 36-41.
6 Algazin S.D. Chislennie algoritmi klassicheskoi mate-maticheskoi fiziki [Numerical algorithms of classical mathematical physics]. Moscow, Dialod-MIFI, 2010, 240 p. (in Russian)
7 Mu L., Wang J., Ye X. Effective implementation of the weak Galerkin finite element methods for the biharmonic equation. Computers & Mathematics with Applications. 2017. vol. 74. no. 6. pp. 1215-1222.
8 Doss L. J. T., Kousalya N. Finite Pointset Method for biharmonic equations. Computers & Mathematics with Applications. 2018. vol. 75. no. 10. pp. 3756-3785.
9 Doss L. J. T., Kousalya N., Sundar S. A Finite Pointset Method for Biharmonic Equation Based on Mixed Formulation. International Journal of Computational Methods. 2017. pp. 1850068.
10 Ryzhskih V.I., Slusarev M.I., Popov M.I. Numerical integration of a biharmonic equation in square area. Vestnik Sankt-Peter-burgskogo universiteta [Bulletin of the Saint-Petersburg university], 2013, no. 10, vol. 1, pp. 52-62. (in Russian)
11 Popov M.I., Soboleva E.A. The approximate analytical solution of the internal problem of conductive and laminar free convection. Vestnik Voronezhskogo Universiteta Ingenernih Tehnologyi [Proceedings of the Voronezh State University of Engineering Technologies], 2016 no. 4, pp. 78-84 (in Russian)
СВЕДЕНИЯ ОБ АВТОРАХ INFORMATION ABOUT AUTHORS
Михаил И. Попов к.ф.-м.н., старший преподаватель, ка- Mikhail I. Popov Cand. Sci. (Phys.-Math.), senior lec-
федра высшей математики и информационных техноло- turer, higher mathematics and information technology de-гий, Воронежский государственный университет инже- partment, Voronezh state university of engineering technol-
нерных технологий, пр-т Революции, 19, г. Воронеж, ogies, Revolution Av., 19 Voronezh, 394036,
394036, Россия,mihail_semilov@mail.ru Russia,mihail_semilov@mail.ru
КРИТЕРИЙ АВТОРСТВА CONTRIBUTION
Михаил И. Попов Полностью подготовил рукопись и несет Mikhail I. Popov Completely prepared the manuscript and is
ответственность за плагиат responsible for plagiarism
КОНФЛИКТ ИНТЕРЕСОВ CONFLICT OF INTEREST
Автор заявляет об отсутствии конфликта интересов. The author declares no conflict of interest.
ПОСТУПИЛА 10.03.2018 RECEIVED 3.10.2018
ПРИНЯТА В ПЕЧАТЬ 02.04.2018 ACCEPTED 4.2.2018
