On the accuracy of difference scheme for Navier-Stokes equations

Full Text

1◦ . Любое численно полученное решение, даже точное аналитическое решение уравнения в частных производных, зависит от ошибок округления, связанных с конечным числом знаков, используемых при арифметических операциях на ЭВМ. Погрешность может оказать существенное влияние на решение конечно-разностных уравнений, так как получение этого решения обычно связано с выполнением большого числа однотипных арифметических операций [1]. В ряде случаев погрешность округления пропорциональна числу узлов разностной сетки, поэтому измельчение сетки, снижая погрешность аппроксимации, может увеличивать погрешность округления. Следовательно, погрешность полученного на ЭВМ решения уравнения в частных производных равна сумме погрешностей аппроксимации и округления. Точность численного решения уравнения в частных производных определяется погрешностью аппроксимации не только самого уравнения, но и граничных условий. Согласованной будет разностная схема, аппроксимирующая уравнеISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1233 © 2014 Самарский государственный технический университет. Образец цитирования: Н. И. С и д н я е в, Н. М. Г о р д е е в а, “О точности разностной схемы для уравнений Навье-Стокса” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2014. № 1 (34). С. 156-167. 156 О точности разностной схемы для уравнений Навье-Стокса ние в частных производных [2]. Условием согласованности разностной схемы является стремление к нулю погрешности аппроксимации при измельчении сетки. Это условие, безусловно, выполняется, если погрешность аппроксимации убывает при измельчении сетки. Для анализа устойчивости конечноразностной схемы применяют, например, метод Фурье [1, 2]. В этом методе решение модельного уравнения представляется рядом Фурье с конечным числом членов и устойчивость (или неустойчивость) определяется тем, что каждое отдельное колебание затухает (или нарастает). Погрешность округления не будет возрастать на каждом шаге по времени, если коэффициент перехода не превосходит единицы. Выражение для коэффициента перехода можно получить при помощи подстановки разложения погрешности в ряд Фурье в конечно-разностное уравнение. В данной статье проводится исследование по времени разностных схем, порядок точности которых может быть произвольно высоким. В работах [1, 2] представлены псевдоспектральные схемы по времени для решения уравнений Навье-Стокса, где использована техника разложений в ряды для поиска особенностей решений уравнений Эйлера [3]. Эти методы обобщены в данной статье на случай схем любого порядка и применения к решению уравнения Бюргерса и уравнений Навье-Стокса [4] в случае несжимаемой жидкости. Предварительно рассмотрим приложение излагаемого метода к уравнению Бюргерса [5, 6], которое является частным случаем уравнений Навье-Стокса в одномерном случае [7, 8]. 2◦ . Рассмотрим случай, когда задана скорость течения жидкости u и её кинематическая вязкость υ . Уравнение Бюргерса [9-13] в общем виде запи˜ сывается так: ∂u ∂u ∂2u +u = υ 2. ˜ (1) ∂t ∂x ∂x Будем искать решение u(x, t) уравнения (1) с начальным условием u(x, 0) = ϕ(x), x ∈ [0; 1], которое является периодическим с периодом 1 по переменной x в любой момент времени t [12, 13]: u(x, t) = u(x + 1, t). (2) Для решения этой задачи посредством разностной схемы третьего порядка по времени применим следующий метод: если известно решение u(x, t) в точке (x, t), то для получения решения в точке (x, t + ∆t) используется аппроксимация третьего порядка, определяемая отрезком ряда Тейлора: u(x, t + ∆t) = u(x, t) + ∂2u (∆t)2 ∂ 3 u (∆t)3 ∂u (x, t)∆t + 2 (x, t) + 3 (x, t) . ∂t ∂t 2 ∂t 6 (3) Частные производные, входящие в формулу (3), вычисляются дифференцированием по времени уравнения (1): ∂u ∂u ∂2u = -u + v 2; ∂t ∂x ∂x 157 Н. И. С и д н я е в, Н. М. Г о р д е е в а ∂2u ∂ = -u 2 ∂t ∂x ∂ ∂3u = -u 3 ∂t ∂x ∂u ∂u ∂u ∂ 2 ∂u - +v 2 ; ∂t ∂t ∂x ∂x ∂t ∂2u ∂u ∂ ∂u ∂ 2 u ∂u ∂2 ∂2u -2 - 2 +υ 2 ˜ . ∂t2 ∂t ∂x ∂t ∂t ∂x ∂x ∂t2 В конце концов, u(x, t + ∆t) можно выразить как функцию от ∆t, u(x, t), (∂u/∂t)(x, t), (∂ 2 u/∂t2 )(x, t), (∂ 3 u/∂t3 )(x, t). Эти производные вычисляются псевдоспектральным методом. Далее изложим обобщение метода [5, 15, 16] на случай схемы произвольно высокого порядка по времени и представим применяемую пространственную схему. Пусть u(x, t) известно при любом x. Тогда аппроксимация порядка n величины u(x, t + ∆t) даётся отрезком ряда Тейлора: u(x, t + ∆t) = u(x, t) + ∂u ∂nu (∆t)n (x, t)∆t + . . . + n (x, t) . ∂t ∂t n! Любую производную (∂ l u/∂tl ) можно получить, дифференцируя (l - 1) раз по времени исходное уравнение (1): ∂lu ∂ l-1 ∂ 2 u ∂ l-1 ∂u u + υ l-1 ˜ . = - l-1 u ∂x ∂x2 ∂tl ∂t ∂t (4) Если предположить решение u(x, t) достаточно регулярным, то можно поменять порядок дифференцирования по x и t. Используя формулу Лейбница, формулу (4) можно записать в виде ∂lu =- ∂tl l-1 i Cl-1 i=0 ∂ i u ∂ ∂ l-1-i u ∂ 2 ∂ l-1 u +υ 2 ˜ , ∂ti ∂x ∂tl-1-i ∂x ∂tl-1 (5) откуда u(x, t + ∆t) = u(x, t)+ n l-1 i Cl-1 - + l=1 i=0 ∂ i u ∂ ∂ l-1-i u ∂ 2 ∂ l-1 u +υ 2 ˜ ∂ti ∂x ∂tl-1-i ∂x ∂tl-1 (∆t)l . (6) l! Искомое решение u(x, t) удовлетворяет (2), что позволяет в любой момент t разложить его в ряд Фурье: +∞ ak (t)ej2πkx , u(x, t) = k=-∞ где 1 j 2 = -1, ak (t) = u(x, t)e-j2πkx dx. 0 Аппроксимация функции u(x, t) определяется методом дискретизации, или псевдоспектральным методом. Точки дискретизации xi определяются 158 О точности разностной схемы для уравнений Навье-Стокса следующим образом: xi = i/N , 0 i < N , где N - число точек дискретизации. Таким образом, используется следующая аппроксимация u∗ (x, t) функции u(x, t): N/2-1 ∗ ak (t)ej2πkx , u (x, t) = k=-N/2 где коэффициенты ak (t) - решения системы уравнений N/2-1 ak (t)ej2πkxi = u(xi , t), 0 i < N, k=-N/2 так как для точек дискретизации выполняется соотношение u∗ (xi , t) = u(xi , t), 0 i < N. Переход от {u(xi , t), 0 i < N } к {ak (t), -N/2 k < N/2 - 1} осуществляется с помощью быстрого преобразования Фурье. Все производные по x вычисляются в пространстве коэффициентов Фурье. Все произведения функций вычисляются в физическом пространстве x. Рассмотрим способ вычисления ∂u/∂x в точке дискретизации xi . Формальная производная функции u∗ имеет вид ∂(u∗ ) (xi , t) = ∂x N/2-1 j2πkak (t)ej2πkxi . k=-N/2 Поскольку в точках дискретизации ∂u ∂x ∗ = ∂(u∗ ) ∂x по определению, то ∂u ∂x ∗ ∂(u∗ ) ∂u (xi , t) = (xi , t) = (xi , t) = ∂x ∂x N/2-1 j2πkak (t)ej2πkxi . k=-N/2 Формула (5) позволяет разработать рекуррентный по l алгоритм вычисления производных ∂ l u/∂tl . Так, например, на l-том шаге известны ∂ i u/∂ti и (∂ i u)/(∂x∂ti-1 ), 1 i l-1, поэтому в правой части формулы (5) неизвестны только две производные: ∂lu ∂x∂tl-1 и ∂ l+1 u . ∂x2 ∂tl-1 Метод [5, 16-18] позволяет получить эти производные посредством двух последовательных дифференцирований по x величины ∂ l-1 u/∂tl-1 , получаемой на (l - 1)-м шаге. После этого величина ∂ l u/∂tl вычисляется как сумма произведений известных величин. Рекуррентный алгоритм позволяет получить всю совокупность производных ∂ l u/∂tl при l, меньших или равных n. После этого по формуле (6) получаем u(x, t + ∆t) и т.д. 159 Н. И. С и д н я е в, Н. М. Г о р д е е в а Описанный процесс позволяет вычислить u(x, t + ∆t) с помощью разностной схемы порядка n, при этом требуется n прямых и 2n обратных преобразований Фурье и используется 2n + 1 массив длины N . Исследование устойчивости осуществлялось для уравнения Бюргерса, линеаризованного в окрестности постоянной скорости U : ∂u ∂u ∂2u +U - υ 2 = 0. ˜ ∂t ∂x ∂x Пусть u(x, t) определена следующим образом: N/2-1 ak (t)ej2πkx . u(x, t) = k=-N/2 Сначала вычисляем производные по времени ∂ l u/∂tl методом [5, 12, 13, 19, 20]: ∂lu (x, t) = ∂tl N/2-1 l ak (t) -j2πkU - υ 4π 2 k 2 ej2πkx = ˜ k=-N/2 N/2-1 = n l -j2πkU - υ 4π 2 k 2 ˜ ak (t) l=0 k=-N/2 (∆t)l j2πkx e . l! Таким образом, n -j2πkU ∆t - υ 4π 2 k 2 ∆t ˜ ak (t + ∆t) = ak (t) l=0 l 1 . l! Разностная схема будет устойчивой в смысле Фурье, если ak (t + ∆t) ak (t) 1 ∀k, - N 2 N - 1. 2 k Положим a = υ 4π 2 k 2 ∆t, b = 2πkU ∆t, тогда схема будет устойчивой, если ˜ n (-jb - a)l gn (a, b) = l=0 1 l! 1 ∀k, - N 2 k N - 1. 2 Заметим, что gn (a, b) - отрезок длины n ряда Тейлора функции e-jb-a и |e-jb-a | = e-a 1, поскольку в данном случае a 0. Более того, при заданных a и b величина gn (a, b) - e-a стремится к нулю при n → ∞. 3◦ . Перейдём к трёхмерной модели. Уравнения Навье-Стокса, определяющие поведение несжимаемой жидкости, в безразмерных переменных записываются следующим образом: ∂u ∂u ∂u ∂u ∂p +u +v +w - υ ∆u + ˜ = 0, ∂t ∂x ∂y ∂z ∂x 160 (7) О точности разностной схемы для уравнений Навье-Стокса ∂v ∂v ∂v ∂v ∂p +u +v +w - υ ∆v + ˜ = 0, ∂t ∂x ∂y ∂z ∂y ∂w ∂w ∂w ∂w ∂p +u +v +w - υ ∆w + ˜ = 0, ∂t ∂x ∂y ∂z ∂z ∂u ∂v ∂w + + = 0. ∂x ∂y ∂z (8) (9) (10) Здесь p - давление; u, v, w - компоненты скорости U ; υ = Re-1 , где Re - ˜ число Рейнольдса. Начальное условие U (x, y, z,0) задано на [0, 1]3 и предполагается периодическим с периодом 1 по переменным x, y, z. Здесь, как и при решении уравнения Бюргерса, аппроксимация величины U (x, t + ∆t), где x = (x, y, z), порядка n записывается в виде отрезка ряда Тейлора U (x, t + ∆t) = U (x, t) + ∂nU (∆t)n ∂U (x, t) ∆t + . . . + n (x, t) . ∂t ∂t n! Основное отличие от предыдущей задачи состоит в том, что уравнения (7)-(9) связаны между собой посредством входящего в них давления и уравнения (10). Предполагая известными величины ∂ i U /∂ti и ∂ i U /∂x∂ti-1 , ∂ i U /∂y∂ti-1 , i U /∂z∂ti-1 , 1 ∂ i l - 1, рассмотрим вычисление ∂ l U /∂tl . Продифференцируем систему уравнений (7)-(10) (l - 1) раз по времени: ∂lu ∂ l-1 ∂ l-1 ∂p = 0, + A - υ l-1 (∆u) + l-1 ˜ ∂x ∂tl ∂t ∂t ∂lv ∂ l-1 ∂ l-1 ∂p = 0, + B - υ l-1 (∆v) + l-1 ˜ ∂y ∂tl ∂t ∂t ∂lw ∂ l-1 ∂ l-1 ∂p = 0, + C - υ l-1 (∆w) + l-1 ˜ ∂z ∂tl ∂t ∂t ∂ l-1 ∂u ∂ l-1 ∂v ∂ l-1 ∂w + l-1 + l-1 = 0, ∂y ∂z ∂tl-1 ∂x ∂t ∂t (11) ∂u ∂u ∂u ∂ l-1 u +v +w , ∂y ∂z ∂tl-1 ∂x ∂ l-1 ∂v ∂v ∂v B = l-1 u +v +w , ∂x ∂y ∂z ∂t ∂ l-1 ∂w ∂w ∂w C = l-1 u +v +w . ∂x ∂y ∂z ∂t (12) где A= Пользуясь формулой Лейбница и предполагая функции u, v, w достаточно регулярными, чтобы поменять порядок дифференцирования по пространству 161 Н. И. С и д н я е в, Н. М. Г о р д е е в а и времени, систему уравнений (12) можно записать так: l-1 i Cl-1 ∂ i v ∂ l-i u ∂ i w ∂ l-i u ∂ i u ∂ l-i u , + i + i ∂ti ∂x∂tl-1-i ∂t ∂y∂tl-1-i ∂t ∂z∂tl-1-i i Cl-1 ∂ i u ∂ l-i v ∂ i v ∂ l-i v ∂ i w ∂ l-i v + i + i , ∂ti ∂x∂tl-1-i ∂t ∂y∂tl-1-i ∂t ∂z∂tl-1-i i Cl-1 A= ∂ i v ∂ l-i w ∂ i w ∂ l-i w ∂ i u ∂ l-i w + i + i . ∂ti ∂x∂tl-1-i ∂t ∂y∂tl-1-i ∂t ∂z∂tl-1-i i=0 l-1 B= i=0 l-1 C= i=0 Векторная запись системы (11) имеет вид ∂lU ∂ l-1 U ∂ l-1 p + Φ - υ ∆ l-1 + grad l-1 = 0, ˜ ∂tl ∂t ∂t div ∂ l-1 U = 0, ∂tl-1 (13) где Φ = (A, B, C). Применяя к уравнению (13) операцию дивергенции, получаем уравнение Пуассона для ∂ l-1 p/∂tl-1 : ∆ ∂ l-1 p ∂lU ∂ l-1 U = - div Φ - div l + υ ∆div l-1 . ˜ ∂tl-1 ∂t ∂t (14) Дифференцируем второе уравнение (13) по t, получаем div ∂lU = 0. ∂tl (15) Из уравнений (13) и (15) следует, что уравнение (14) можно записать в виде ∆ ∂ l-1 p = - div Φ. ∂tl-1 (16) Уравнение Пуассона (16) решается спектральным методом [18-20]. После этого достаточно найти grad(∂ l-1 p/∂tl-1 ), чтобы из первого уравнения (13) получить ∂ l U /∂tl . Следует ещё раз отметить, что решение рассматриваемой задачи периодично в направлениях x, y, z с периодом 1. Таким образом, в любой момент времени t можно разложить U (x, t) и p(x, t) в ряды Фурье: +∞ +∞ +∞ U (x, t) = ak(t)ej2πkx , kx =-∞ ky =-∞ kz =-∞ +∞ +∞ +∞ bk (t)ej2πkx , p(x, t) = kx =-∞ ky =-∞ kz =-∞ k = (kx , ky , kz ), 162 x = (x, y, z), j 2 = -1, kx = kx x + ky y + kz z. О точности разностной схемы для уравнений Навье-Стокса Пусть Ω - множество точек дискретизации: Ωx = {i/Nx , 0 Ω = Ωx × Ωy × Ωz , i < Nx } , Ωy = {i/Ny , 0 i < Ny } , Ωz = {i/Nz , 0 i < Nz } . Следует отметить, что в процессе решения использовалось допущение, что уравнение divU = 0 остаётся справедливым на Ω в каждый момент времени t. В рассмотренных ниже случаях при небольшом числе шагов по времени величина divU / U остаётся меньшей 10-8 . С целью исследования устойчивости разностных схем для уравнений Навье-Стокса рассматривались уравнения, линеаризованные в окрестности постоянной скорости U0 = (u0 , v0 , w0 ). Предполагалось, что U (x, t) и p(x, t) определены следующим образом: Nx /2-1 Ny /2-1 Nz /2-1 U (x, t) = ˆ Uk (t)ej2πkx , kx =-Nx /2 ky =-Ny /2 kz =-Nz /2 Nx /2-1 Ny /2-1 Nz /2-1 pk (t)ej2πkx , ˆ p(x, t) = kx =-Nx /2 ky =-Ny /2 kz =-Nz /2 ˆ Uk (t) = uk (t), vk (t), wk (t) . ˆ ˆ ˆ Чтобы U и p были решениями линеаризованных уравнений, необходимо, чтобы коэффициенты Uk (t) и pk (t) удовлетворяли соотношениям ˆ ∂ uk ˆ + j2πkx pk + uk j2π(kx u0 + ky v0 + kz w0 ) + υ 4π 2 k 2 = 0, ˆ ˆ ˜ ∂t ∂ˆk v + j2πky pk + vk j2π(kx u0 + ky v0 + kz w0 ) + υ 4π 2 k 2 = 0, ˆ ˆ ˜ ∂t ∂ wk ˆ ˜ + j2πkz pk + wk j2π(kx u0 + ky v0 + kz w0 ) + υ 4π 2 k 2 = 0, ˆ ˆ ∂t kx uk + ky vk + kz wk = 0. ˆ ˆ ˆ (17) Из уравнений (17) следует, что 2πk 2 pk = j kx ˆ ∂ uk ˆ ∂ˆ v ∂w ˆ + ky k + kz k ∂t ∂t ∂t . Но метод, изложенный в [5, 6], предполагает, что ∂ U /∂t удовлетворяет уравнению несжимаемости: div ∂U = 0, ∂t kx ∂ uk ˆ ∂ˆ v ∂w ˆ + ky k + kz k = 0. ∂t ∂t ∂t Давление p удовлетворяет уравнению ∆p = 0 при условиях периодичности. Получаем, что давление постоянно и не входит в дальнейшие вычисления. Поэтому уравнения (17) разделяются, и исследование устойчивости 163 Н. И. С и д н я е в, Н. М. Г о р д е е в а сводится к исследованию устойчивости следующего уравнения: ∂ uk ˆ + uk j2π(k U0 ) + υ 4π 2 k 2 = 0. ˆ ˜ ∂t Согласно разностной схеме, описанной в [2, 7, 8], u(x, t + ∆t) определяется коэффициентами uk (t + ∆t) отрезка ряда Фурье, удовлетворяющими соотноˆ шению n l 1 . uk (t + ∆t) = uk (t) ˆ ˆ -j2π(k U0 )∆t - υ 4π 2 k 2 ∆t ˜ l! l=1 Поэтому вновь воспользуемся функциями gn , определёнными в [12-14]. Схема устойчива, если -Nx /2 gn υ 4π 2 k 2 ∆t, 2π(kx u0 + ky v0 + kz w0 )∆t ˜ 1 ∀k, kx Nx /2 - 1, -Ny /2 ky Ny /2 - 1, -Nz /2 kz Nz /2 - 1. 4◦ . Рассматривалось решение уравнения (1) с начальным условием u(x, 0) = - sin(2πx) (18) и условием (2) при υ = (2π)-1 , которое имеет следующий аналитический вид: ˜ u(x, t) = 4 a0 + +∞ -2πp2 t sin(2πpx) p=1 pap e , 2 +∞ ap e-2πp2 t cos(2πpx) p=1 ap = (-1)p Ip (1/2), где Ip - модифицированная функция Бесселя первого рода. Значения отклонений расчетного ucalc от точного uexact вычислялись в норме l2 для N -мерного пространства в момент времени t = 0.02 для различных значений порядка разностной схемы n: 1 ∆= N 1/2 N -1 uexact (i/N, 0.02) - ucalc (i/N, 0.02) 2 . (19) i=0 Результаты вычисления ∆ для решения, полученного по разностной схеме, представлены на рисунке сплошными линиями. На этом же рисунке приведены результаты (штриховые линии), полученные в программной среде Maple 6. Отсутствие некоторых значений ∆, вычисленных в Maple 6, объясняется тем, что либо вычисления длятся очень долго, либо данный алгоритм в Maple 6 расходится. В случае N = 8 ошибка ∆ убывает при росте порядка схемы от 1 до 4 и остается постоянной для порядков, больших 4. Таким образом, при N = 8 нельзя надеяться получить хорошую точность схем. В случае N = 16 то же явление имеет место при порядке схем более семи. При N = 32 и N = 64 имеем очень хорошую точность. Проведенные вычисления показали, что в рассмотренных случаях точность вычислений заметно возрастала при возрастании порядка разностной 164 О точности разностной схемы для уравнений Навье-Стокса Кривые отклонения (19) для уравнения (1) с начальным условием (18) и условием (2) при t = 0.02: 1 - N = 8, ∆t = 2 · 10-2 ; 2 - N = 16, ∆t = 5 · 10-3 ; 3 - N = 32, ∆t = 1,25 · 10-3 ; 4 - N = 64, ∆t = 3,125 · 10-4 ; сплошные линии - вычисление по разностной схеме; штриховые линии - вычисленные в программной среде Maple 6 Deflection curves (19) for equation (1) with initial condition (18) and condition (2) when t = 0.02 for N = 8, ∆t = 2 · 10-2 (curves number 1); N = = 16, ∆t = 5 · 10-3 (curves number 2); N = 32, ∆t = 1,25 · 10-3 (curves number 3); N = 64, ∆t = 3,125 · 10-4 (curve number 4). Solid curves are obtained by the difference calculating scheme. Dashed curves are obtained in Maple 6. схемы по времени. Однако при этом не учитывались издержки, связанные со временем использования центрального процессора для расчетов. Показано, что алгоритмы быстрого преобразования Фурье позволяют создавать спектральные и псевдоспектральные методы аппроксимации в пространстве. Полученные методы могут использоваться в конечно-разностных схемах, которые имеют высокий порядок аппроксимации в пространстве, но сравнительно низкий порядок точности по времени. Показано, что решение уравнений в частных производных псевдоспектральным методом позволяет получить весьма высокую точность при вычислении пространственных производных. Результаты работы показывают, что конечноразностная схема, используемая для вычисления производных по времени, является основным источником отклонений приближенного решения от точного решения.

About the authors

Nikolay I Sidnyaev

N. E. Bauman Moscow State Technical University

Email: sydnyaev@bmstu.ru
(Dr. Techn. Sci.), Head of Department, Dept. of Higher Mathematics 5/1, 2-ya Baumanskaya st., Moscow, 105005, Russian Federation

Nadezhda M Gordeeva

N. E. Bauman Moscow State Technical University

Email: nmgordeeva@bmstu.ru
Assistant, Dept. of Higher Mathematics 5/1, 2-ya Baumanskaya st., Moscow, 105005, Russian Federation

References

Supplementary files