Численное решение линеаризованной системы Осколкова Текст научной статьи по специальности «Математика»



Серия «Математика»

2015. Т. 11. С. 23—34

Онлайн-доступ к журналу: http://isu.ru/izvestia

УДК 517.9

Численное решение линеаризованной системы Осколкова *

П. Н. Давыдов1, М. В. Плеханова2'3

1 Федеральный центр сердечно-сосудистой хирургии, 2Южно-Уральский государственный университет, 3 Челябинский государственный университет

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

Ключевые слова: жидкость Кельвина - Фойгта, система уравнений Осколкова, численное решение начально-краевой задачи, устойчивость разностной схемы, сходимость разностной схемы.

Рассмотрим начально-краевую задачу для системы Осколкова динамики вязкоупругой жидкости Кельвина - Фойгта в ограниченной

области О с Кп с границей дО класса С^

(1 - хД)™г = VДад - (IV ■ У^ -Ур, (х,г) е О х [0,Т], (0.1)

Здесь параметр х е К, как правило, характеризует упругие свойства жидкости, а параметр V е К — её вязкие свойства. Вектор-функции и> = (■ш1,,ш2,...,,шп) (вектор скорости жидкости), Ур = г = (г1,г2,...,гп) (градиент давления) неизвестны. Системы уравнений такого вида

* Работа поддержана РФФИ (грант 14-01-31125 мол_а)

Введение

У--ш = 0, (х,г) е О х [0,т], Цх,*)=0, (х,г) е дО х [0,Т], ■ш(х, 0) = ■ш0(х), х е О.

(0.2) (0.3)

(0.4)

встречаются в динамике неньютоновских жидкостей [5; 2] и относятся к системам, не разрешимым относительно производной по времени. Они не укладываются в рамки классической математической физики и требуют отдельных методов исследования, качественного и численного [1; 4].

Разрешимость задачи (0.1)-(0.4) в различных классах решений исследовалась в работах [5; 2; 7]. В данной работе предложена разностная схема для численного решения линеаризованной задачи, доказана ее устойчивость и сходимость решения разностной системы уравнений к решению исходной задачи. При этом параметры задачи выбирались в соответствии с условиями, гарантирующими ее однозначную разрешимость согласно работе [7].

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

Численная схема строилась с помощью замены производных в системе на их разностные аналоги. Однако вырожденность системы потребовала нетривиального вывода разностных соотношений для нахождения давления, производная по времени от которого отсутствует в системе. Аппроксимация доказывалась с помощью разложений по формулам Тейлора разностных функций, а устойчивость схемы — на основе метода разделения переменных [3].

Выведенная схема численного решения была реализована в программном продукте [6], позволяющем проводить численный поиск решений задачи, изучать решения при различных параметрах, осуществлять численный эксперимент. Программная часть работы реализована на языке программирования С# в виде обьектно-модульного программного комплекса. Использованы свободные библиотеки для построения 2Б, 3Б графиков для работы с матрицами. Разработанная разностная схема и методы исследования могут в дальнейшем стать отправной точкой при исследовании близких по структуре математических моделей.

1. Линеаризованная модель

Рассмотрим начально-краевую задачу для линеаризованной системы Осколкова в окресности стационарного решения «(х),и(х)):

(1 -хД)«г = vДw - (ги -V)« - (г -V)«; -Ур, (x,í) е П х [0,Т], (1.1)

V-« = 0, (х,г) е П х [0,Т], (1.2)

г(х,Ь)=0, (х,г) е дП х [0,Т], (1.3)

г(х, 0) = г0(х), х е П. (1.4)

Здесь П с Мга — ограниченная область с границей с дП класса С^.

Выберем п = 2, г = (и, V), Ур = (рх,ру), П = (0,п) х (0,п). Тогда система уравнений (1.1)—(1.4) принимает вид

иг - хигхх - хиуу = Vихх + ииуу--иих - ■ииу - иих - vuy - рх, (х, у, ¿) е П х [0, Т], (1.5)

VI - XVtxx - XVtyy = VVxx + VVyy -

-тх - ту - ииих - vUy - ру, (х, у, ¿) е П х [0, Т], (1.6)

их + Vy = 0, (х,у,г) е П х [0,Т], (1.7)

и(х,¿) = 0, v(x,t) = 0, (х,г) е дП х [0,Т], (1.8)

и(х, 0) = и0(х), v(x, 0) = vo(x), (х,у) е П. (1.9)

Введем равномерную сетку на области П с шагом Н = п/К, определив тем самым точки х1 = гН, у^ = ,]Н, г = 0,...,К, ] = 0,...,К. Определим шаг по временни т > 0, получив точки Ьп = пт разбиения отрезка [0,Т], п = 0,...,Ж. Приближенные значения функций и^,р в узлах с координатами (х^,у^,Ьп) будем обозначать через и'^, иП^,р1]. Построим разностную схему, используя левые аппроксимации производных для функций и^ в уравнении (1.7) и для функции р в уравнениях (1.5), (1.6), правые аппроксимации для всех остальных производных в уравнениях (1.5), (1.6).

Для краткости в формулах возьмем и = 0, и = 0 и введем обозначения

= «П+и + г'П-и + +1 + «П— - 4«?-,

у,ру) = + ^А^ - ^(рЬ+1 ~ РЬ)-

Тогда разностная схема примет вид

-^КГ + = (1-Ю)

Т п п п п п 1 п п п п

ут-Х, №

,п I Л„,п I „п

\ __I_ „«■ I_ „«■ I_ „«■ _ _ „«■ _ „«■ _ „«■ _

-иг-1,]-1 + 4иг-1+ у1+1,3 + ^г—1'] + '^п,]+1 +

+<— - 4'п] - ''п+1]-1 - 'п-1]-1 - уп - уп-2 + -1] (1.12)

с начальными условиями

и°ц = и0'П], = УО'П]

и граничными условиями

уп.м = 0,

и0,3 = иМ,] = ип,0 = = '0,3 = 'N,3 = иг,0 =

Ро,3 = Рм,з = Рпо = Рпм = 0, 1,3 =0,1,...,К, п = 0,1,...,Ы.

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

2. Порядок аппроксимации численного метода

Выразим невязку в узлах сетки (1.10)-(1.12) для линеаризованной модели Осколкова:

~^2^Уг+1,3 ^ иг~1,3 ^ УгЛ+1 ^ ~ ^УгЛ ^

+ иг-1,3 + Уг,3+1 + Уг,3~1 ~ ^иг,з)~

~(уТ?+1,з + ,3 + уХз+1 + уХз-1 ~~ 4-и™,-), (2.1)

П = \(Р1з+1-Р1з) + \(у^1-<з)-

~^УШз + + + 1 - +

+ ^2 (""¿+1,3 + ,3 + У^3 +1 + УХз~ 1 ~~

Известия Иркутского государственного университета. 2015. Т. 12. Серия «Математика». С. 23-34

ЧИСЛЕННОЕ РЕШЕНИЕ ЛИНЕАРИЗОВАННОЙ СИСТЕМЫ ОСКОЛКОВА 27 + + <7+1 + " (2-2)

т п п п п п 1 п п п п

vт - X

Н3

Гч<п I „п I „п , п д п п п _ п

№+1 ,з + иг-1 ,з + и,г,з+1 + а1,з-1 Ааг,з а1,з и-2,з а—1 ,з+1

-ип-1 ,з-1 + 4ип-1 ,з + ^+1 ,з + vn-1 ,з + ^,3+1 +

+^^-1 - АУпз - ^++1,3-1 - - VIпз - VIпз-2 + -1]. (2.3)

С помощью формулы Тейлора приведем невязку (2.1) к виду

С = Рх + ^Рхх + о{К) +Щ + + о(т) —

-Х(игхх + игуу + тиггхх + тиЫуу) + о(т + Н) - V(ихх + иуу) + о(Н)

Нт

= ^Рхх + ~ ХТ(иих.

Аналогично из (2.2) получим

Нт

= -^Рхх + 2 м« ~ Хт(Щгхх + 1Чгуу) + о(т + й) = О(т) + 0(/г).

Нт Я = Ру + т;Руу + ^ + 2

-Х^гхх + ^^гуу + тvttxx + ^иуу) - V^хх + Vyy) + о(т + Н) = Нт

= 2Руу + ~ Хт(Уахх + ^УУ) + + = °(г) +

Рассмотрим невязку для уравнения (2.3) 61 = т 7о(Р хххх + Руууу) + О (Л. ) ) + (ихх + 1*^) +

/ Н \

ит) I ~{иХххх 11ххуу УХХуу Ууууу) о{К) ) о{К) —

= О(Н) + О(тН).

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

Др - ° иххх + ихуу + "^хху + "Уууу - Д(их + 'Уу) - 0.

Таким образом, если точное решение и^,р задачи (1.5)-(1.9) при и = и = 0 четырежды непрерывно дифференцируемо, то разностная схема (1.10)-(1.12) аппроксимирует задачу с первым порядком по т и Н.

3. Устойчивость разностной схемы

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

К = Е р>щх^, < = Е <еЩхгегау, Р = Е егау1.

Яув

Подстановка в разностную схему после несложных преобразований приводит к соотношениям

(Pqs Pqs)

Т

"¿2 - 2p™+1 + pngtle-%qh ~ (pnqse^h - 2pnqs + р™ е"^))-

" + Pnqfl?~lsh ~ (Pnqaeiah ~ 2 Pnqs + Pqse~ish)) =

= Лр, (pqse<1 — 2 Pqs + pqse q ) + (Pqse ~ 2Pqs + Pqse

1 h 1

■I О

qs y qs

/ _n \

"^2 - 2<*+1 + - - 2anqs + onqse-^h))~

X i „.n+1 Ash o.-n+1 i „.n+1„-ish / n Ash n'n , n —ish\\ (aqs e ~ qs + <V в ~ (aqse ~ 2(7qs + aqse )) =

= ~ 2+ °>~гф) + ~ 2+ <se~lsh)-

(3.2)

n iqh n n ish n Pcsse q - P(is + V3 - V _ 0 (3'3)

Формула давления после преобразований дает выражение

Kn eiqh i Kn pish + Kn p-iqh + Kn p-ish _ »Kn _ Kqsp + Kqsp + Kqsp + Kqsp 4Kqs _

h

7„n „n „-iqh* n n —ish\ , = ~ \ Pqs ~ Pqse + ~ ) +

Т

+ p™ e"^ + p" eisfc + p™ e"^ - 4p™ -

_ pn p-2iqh _ pn p-iqhpish _ n p-iqhp-ish + 4pn p-iqh+

qs - qs - qs - qs qs

+^гф + + °>гзН + - 4опяз-

-о°3еПдкв-пз1г - о°3е-Пдке-П8к - о°8 - од°3е-2П8Н + 4о°3е-П8к]. Теперь выражения (3.1), (3.2), (3.3) можно переписать в виде

\ (Р^1 - Я.) - % ("4р™+1 «ш2 ^ + 4р-3 8ш2 ^ ) -

— ¿Vьш 2 Л,2 43 2 ¡г \ ^ ^) ' \ - - 1 ЙП2 ^ + Ф7» 8Ш2 -

Х (-Аа^1 8т2 ^ + 4а" 8т2 =

Тк2 у ~д8 2 , ^д^ 2 ^

-4v п . 2 дН -4v п . 2 вН 1 / п 8Ь п -а1 8Ш2 + -ГТ5-0-" 8Ш2 — - - Сем/г - к?

Н2 д8 2 ' Н2 2 НЧ д8 д8

,п+1 / 1 , 4Х Л!т,2 Ф , ^

д8

= (7 + (*,» | + 1)) - 1к;.„ (3.5)

, п ( . 2 дН . 2 вН

( 8Ш2 у + 81П2 у

= ^ (1 " е"^) + <т£ (1 - е"^) ) +

VТ - х

2 аН 2 вН\

2 +ВШ т)

-4 ( 8Ш2 V + 8Ш2 — ) (Р™ + (7™ )-

-Р°8(1 + е-2пдН + е-ПдНеП8Н + е-пфе-П8Н - 4е-пф)--о°8(1 + е-2П8Н + еПдНе-П8Н + е-ПдНе-П8Н - 4е-П8Н) \ . (3.6)

Выделим главные члены ассимптотического разложения по малому параметру Н выражений по обе части равенства (3.6) и получим

„ п (а2Н2 а4Н4 в2Н2 в4Н4

— Лкп I _--_--1____

дз\ 4 48 4 48

= 7 + + а"8 + Нг)) +

+

ут-Х

кт

21,2 , д4Н4 2 1,2 , в4Н4\ (

п + оп д8 + о д8 -

7а4Н4 в4Н4 д2в2Н4

-р™ -к2{д2 + ¿2) + 1дАкА + 1дв к + —--1---Ь

12

12

-С ( + *2) + гз^ + гд2з^ + + + ^^

12

12

2

Г)

+

VТ - х

з41г2\ 1

12 у ~т

'П /

рд 8 -гд3

2 д8

Ли „2^2

2 д4Н д2в2Н

+опЛ -гв3 - гд2в -

22

в4Н д2в2Н\

+

кп -__

22

т (д2 + в2)\ 12(д2 + в2)

п Л а2Н г ,( -3-2 а4Н а2в2Н\\ +Рд8 [гд + у- + (ут - X) [~гд3 - гдз2 - - ) J +

п Л в2Н . ,3.2 в4Н д2в2Н\\

[г8 + — + т -X) - п в - — - ) ) ,

+ <%„ ( + +

2

кп =

дз т{д2 + 52)

Ша(^т - х)(д2 + в2) -1)) +

+опм^т - х)(д + в2) -1) =

V - х)(д2 + в2) -1,,

(гдрп°3 + гво%). (3.7)

Т(д2 + в2)

Используя (3.7), рассмотрим асимптотическое представление при малом Н соотношения (3.4)

«+11Л , (IИ. _ , _ ^

Р«3 - 1 + К2 V 4 48 4 48 У У

1(л 4(х - vт) (д2Н2 д4Н4 в2Н2 в4Н4\\

= РпЯ8~ ( 1 +

\2 I

48

+

48

+

+

V - х)(д2 + в2) -1, 2

{д2рп°8 + двоп8,

рп+1 = д8

д8

1 + х(д2 + в2)

Т (д2 + в2)

1 , / М 2 . 2\ , V - х)(д2 + в2) - 1 2\ . 1 + (х-"т)(д + +- + -д ] +

и

1

4

4

п

+ 1+ ^.Л-Ч'8^-1«'. (3-8)

1 + х(Ч2 + 82) Ч2 + 82 Аналогично для (3.5) получим

= 'Ь С! + (х _ „^ + ^ + +

д3 1+ Х(Ч2 + 82) V Ч + 82

Положим К = ч2 + в2. Полученные соотношения (3.8), (3.9) можно записать в виде

рп+^\ ( А в\( рп

^п+ч V С Б )\о%а

(3.10)

где

(ит-х)К-1 {ут-Х)К- 1

"Щ+х^Г9 ' ЛГ(1 + 9*'

Характеристическое уравнение матрицы (3.10) имеет вид

(А - А)(Б - А) - ВС = А2 - (А + Б)А + АБ - ВС = 0, нетрудно вычислить выражения

1+ хК 1 + хК

Ав = {1 + хку[1 + 2(х " ит)К + (х " ит)2кЧ

+(1 + (x - »т)кж„г - х)к -1) + (г - л2*2! _ - x) v*2

К2

К2(1 + хК)2 '

К 2(1 + хК )2

Поэтому собственными значениями матрицы (3.10) являются числа

vK

А = 0, \ = А + Б = 1 — ————т.

1 + хК

Сходимость численного метода линеаризованной модели Осколкова гарантируется при выполнении условия |А| < 1, которое в случае х > 0, V > 0 выполняется, когда 1 — > —1, т. е. г < 2(-1г^Х->. Это в

любом случае верно при г < Доказано следующее утверждение.

Теорема 1. Пусть х > 0, V > 07 т < Тогда при достаточно малом Н схема (1.10)-(1.12) устойчива.

Замечание 1. Строгое неравенство |А| < 1 использовано, потому что для вычисления собственных значений А матрицы использовалась асимптотика по малым Н.

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

4. Численный эксперимент

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

1. Задаются границы пространственной переменной, длина шага по времени. Разбивается заданный интервал на части. Для функций u и v задаются начальные и граничные значения, а также параметры v, %.

2. Методом матричной прогонки по известным значениям нулевого временного слоя функций u и v находятся значения на нулевом слое функции p. После подстановки известных значений функций u,v,p на нулевом временном слое, получаются значения на первом временном слое функций u, v. Продолжается поочередное вычисление значений функций u, v, p.

Пусть х = 0.001, v = 0.001, начальные значения u0 = sin2 x sin y cos y, v0 = - sin2 y sin x cos x, T = 0.01, Mi = 50, M2 = 50, h = 0.06.

На рисунке ниже приведен график численного решения на 50-м шаге по времени.

■Р[х,50] -13,6/0,5

У(х.50] -ОАЭ / 0,51

ЧИСЛЕННОЕ РЕШЕНИЕ ЛИНЕАРИЗОВАННОЙ СИСТЕМЫ ОСКОЛКОВА 33 Список литературы

1. Демиденко, Г. В. Уравнения и системы, не разрешенные относительно старшей производной / Г. В. Демиденко, С. В. Успенский. - Новосибирск : Науч. книга, 1998. - 436 с.

2. Звягин, В. Г. Исследование начально-краевых задач для математических моделей движения жидкостей Кельвина - Фойгта / В. Г. Звягин, М. В. Турбин // Современная математика. Фундаментальные направления. - 2009. - Т. 31. -С. 3-144.

3. Калиткин, Н. Н. Численные методы / Н. Н. Калиткин. - СПб. : БХВ-Петербург, 2011. - 592 с.

4. Линейные и нелинейные уравнения соболевского типа / А. Г. Свешников, А. Б. Альшин, М. О. Корпусов, Ю. Д. Плетнер. - М. : Физматлит, 2007. -734 с.

5. Осколков, А. П. Начально-краевые задачи для уравнений движения жидкостей Кельвина - Фойгта и жидкостей Олдройта / А. П. Осколков // Тр. Мат. ин-та АН СССР. - 1988. - Т. 179. - С. 126-164.

6. Свидетельство о государственной регистрации программы для ЭВМ № 2014618458; РОСПАТЕНТ/ П. Н. Давыдов, М. В. Плеханова; заявитель и правообладатель ФГБОУВПО «Челябинский государственный университет». - № 2014616031; заявл. 24.06.2014; гос. регистрация в Реестре программ для ЭВМ 20.08.2014.

7. Федоров В. Е. Полулинейные вырожденные эволюционные уравнения и нелинейные системы гидродинамического типа / В. Е. Федоров, П. Н. Давыдов // Тр. Ин-та математики и механики УрО РАН. - 2013. - Т. 19, № 4. - С. 267-278.

Давыдов Павел Николаевич, Федеральный центр сердечно-сосудистой хирургии, 454003, г. Челябинск, проспект Героя России Родионова Е. Н., 2 тел.:(351)7342795 (e-mail: davydov@csu.ru)

Плеханова Марина Васильевна, кандидат физико-математических наук, доцент, Южно-Уральский государственный университет (НИУ), 454080, Челябинск, пр. Ленина, 76 тел. (351)2679900 Челябинский государственный университет, 454001, Челябинск, ул. Бр. Кашириных, 129 тел. (351)7997235 (e-mail: mariner79@mail.ru)

P. N. Davydov, M. V. Plekhanova

Numerical Solution of the Linearized Oskolkov System

Abstract. Difference scheme is constructed for the numerical solution of the linearized Oskolkov model for Kelvin - Voigt fluid. The approximation of the scheme with the first order, stability and convergence has been proven. Problem-oriented programs complex is worked out for the numerical solution of the corresponding problem. Computing experiment was realized by means of the complex.

Keywords: Kelvin - Voight fluid, Oskolkov system of equations, numerical solution of initial boundary value problem, difference scheme stability, difference scheme convergence.

References

1. Demidenko G.V., S.V. Upsenskii. Partial Differential Equations and Systems not Solvable with Respect to the Highest-Order Derivative. New York, Basel: Marcel Dekker, 2003.

2. Zvyagin, V. G., Turbin M. V. The study of initial-boundary value problems for mathematical models of the motion of Kelvin-Voigt fluids. Journal of Mathematical Sciences, 2010, vol. 168, № 2, pp.157-308.

3. Kalitkin N.N. CHislennye metody [Numerical methods]. SPb., BHV-Peterburg, 2011, 592 p. (in Russian)

4. Sveshnikov A. G., Al'shin A. B., Korpusov M. O., Pletner JU. D. Linejnye i nelinejnye uravnenija sobolevskogo tipa [Linear and Nonlinear Equations of Sobolev Type]. M., Fizmatlit, 2007, 734 p. (in Russian)

5. Oskolkov A.P. Initial-boundary value problems for equations of motion of Kelvin - Voight fluids and Oldroyd fluids. Proceedings of the S'teklov Institute of Mathematics, 1989, vol. 179, pp. 137-182.

6. Svidetel'stvo o gosudarstvennoj registracii programmy dlja [Certificate of state registration of the computer number] JEVM № 2014618458; ROS-PATENT/ Davydov P.N., Plehanova M.V.; zajavitel' i pravoobladatel' FGBOUVPO «Cheljabinskij gosudarstvennyj universitet». [the applicant and the right holder FGBOUVPO «Chelyabinsk State University».] № 2014616031; zajavlenie [Statement] 24.06.2014; gosudarstvennaja registracija v Reestre programm dlja JEVM [state registration in the Register of Computer Programs] 20.08.2014. (in Russian)

7. Fedorov V.E., Davydov P.N. Polulinejnye vyrozhdennye jevoljucionnye uravnenija i nelinejnye sistemy gidrodinamicheskogo tipa [Degenerate semilinear evolution equations and nonlinear systems of hydrodynamic type]. Trudy instituta matematiki i mehaniki UrO RAN. [Proceedings of the Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences], 2013, vol. 19, no 4, pp. 267-278. (in Russian)

Davydov Pavel Nikolayevich, Federal Center for Cardiovascular Surgery, 2, Av. Hero of Russia Rodionov E.N., Chelyabinsk, 454003 tel.: (351)2679900 (e-mail: davydov@csu.ru)

Plekhanova Marina Vasilyevna, Candidate of Sciences (Physics and Mathematics), National Research South Ural State University, 76, Av. Lenin, Chelyabinsk, 454080 tel.: (351) 267-99-00 Chelyabinsk State University, 129, Kashirin Brothers st., Chelyabinsk, 454001 tel.: (351) 799-72-35 (e-mail: mariner79@mail.ru)