Методика построения разностных схем для задачи диффузии-конвекции-реакции, учитывающих степень заполненности контрольных ячеек Текст научной статьи по специальности «Математика»

УДК 519.6

А.И. Сухинов, А.Е. Чистяков, Н.А. Фоменко

МЕТОДИКА ПОСТРОЕНИЯ РАЗНОСТНЫХ СХЕМ ДЛЯ ЗАДАЧИ ДИФФУЗИИ-КОНВЕКЦИИ-РЕАКЦИИ, УЧИТЫВАЮЩИХ СТЕПЕНЬ ЗАПОЛНЕННОСТИ КОНТРОЛЬНЫХ ЯЧЕЕК

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

Дискретная модель; погрешность аппроксимации; устойчивость; консервативность; уравнение диффузии-конвекции-реакции.

A.I. Sukhinov, A.E. Chistyakov, N.A. Fomenko

METHOD OF CONSTRUCTION DIFFERENCE SCHEME FOR PROBLEMS OF DIFFUSION-CONVECTION- REACTION, TAKES INTO THE DEGREE FILLING OF THE CONTROL VOLUME

The authors propose a variant of the integro-interpolation method that considers the degree of fullness of the control cells. The presented approximations are proven to be the second order on the spatial variable. The paper contains description of the discrete mathematical model and proof of its stability, based on the mesh maximum principle. The proposed mesh equations are written in the canonical form. We also check conditions that must be fulfilled for the mesh maximum principle was applicable and show conditional stability of the difference equations by initial data, boundary conditions and left part. The operator of the diffusion transfer is shown to be conservative.

Mesh equations; diffusion-convection-reaction equation; approximation; stability; conservativity.

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

Постановка задачи. Рассмотрим задачу транспорта веществ, которая может быть представлена уравнением диффузии-конвекции-реакции

< + uc'x + v< = (К )'x + (К )'y + f (1)

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

сп (х, у, О = апс+Рп, (2)

где ы,у - компоненты вектора скорости; ^ - коэффициент турбулентного обмена; / - функция, описывающая интенсивность и распределение источников.

Построение дискретной модели. Расчетная область вписана в прямоугольник. Для численной реализации дискретной математической модели поставленной задачи вводится равномерная сетка:

^ ={гп = пт, х. = 1НХ, у. = ]Ну; п = 0...Ы,, г = 0... Их, ] = 0...Иу;

М,т = Т, Мхкх = 1х, Муку = 1у},

где т - шаг по времени; кх, ку - шаги по пространству; - верхняя граница по

времени; Их, Иу - границы по пространству.

Для аппроксимации уравнения (1) по временной координате используем схемы с весами [1-2]

c - c

+ iicx+ Vcy ={ßcx )x+(ßcy) + f ,

(3)

т V - 'у

где с = ас + (1 -а)с , <ге [0,1] - вес схемы.

Ячейки представляют собой прямоугольники, они могут быть заполненными, частично заполненными или пустыми.

Рис. 1. Расположение ячейки относительно прилегающих к ней узлов

Центры ячеек и узлы разнесены на кх/2 и Ну12 по координатам х и у соответственно. Обозначим через ог. заполненность ячейки (г, ]). Поле концентрации рассчитывается в вершинах ячейки, как показано на рис. 1. Вершинами ячейки (г, ]) являются узлы (г, ]), (г -1, ]), (г, ] -1), (г -1, ] -1). В окрестности узла (г, ]) лежат ячейки (г, ]), (г +1,.), (г, ] +1), (г +1, ] +1), как показано на рис. 2.

Вводятся коэффициенты я0, д1, я2 , я3, я4 , описывающие заполненность контрольных областей, находящихся в окрестности ячейки [3]. Значение 40 характеризует заполненность области Б0: хе (-1/2,х!+1/2), уе(у.-1/2,у.+1/2), цх - Д:

х е(х , х1 + 1/2 ) у е ( (.-1/2 , у. +1/2 ) 42 - °2: х е (хг-1/2 , хг ), у е (Уj-1/2, у.+1/2 ) 4з - А:

Х е VХ' -1/2' Х' +1/2 )' У е( ' У^+1/2 )' Ча, : хе (-1/2'х{+1/2), у е(уj-1/2'у]). Запол-

ненные части областей Бт будем называть От, где т = 0...4. В соответствии с этим коэффициенты чт можно вычислить по формулам

. . S0 о. .+ о. +, .+ о.+, .+, + о. .+, ч о. +, .+ о.+, .+,

(Чт )„ = ^ ' (Чо )„ = "j ' + 1'j 4' + 1' j + 1 " j + 1 ' (Ч )„ = ' + 1'j 2 ' + 1'j + 1 ' вт

о ■ + о .+, ч о.+, .+, + о. .+, о. . + о.+, .

(Ч2 )и = ' (Чз)'', = ' (Ч4)'', = .

Рис. 2. Расположение узлов относительно ячеек

Проинтегрируем по области О0 уравнение (3) и воспользуемся свойством линейности интеграла' в результате чего получим: •с - с

Ц-ёхёу + Ц исХ йхду + Ц усу дхйу =

Оо Оо Оо

/ ' Ц (/ис'х )х dxdу + Ц (/ис'у) dxdу + Ц fdxdу . (4)

Оо Оо Оо

Вычислим отдельно каждый из полученных интегралов:

Л — (Чо)''. Ц = (Чо)''; ЪЪу. (5)

Оо

Второй интеграл в выражении (4) запишется так: Цисхdxdу = Цисdxdу + Цисхdxdу - (ч1 ) Цисхdxdу + (Ч2 ) j Цисхdxdу .

Оо О1 О2 и1 в2

Вычисляем интегралы по областям Д и Б2:

у,+1/2 х: +1/2 с — с

Ц исХdxdу = | dу | исХdx — имп ] '+1''-— Ну'

п 2

y j+1/2

Ц ucxdxdy = J dy J ucx dx — u

cK j - ci -i,j

i-1/2, j

D2 y j-1/2 xi-1/2 в результате чего получим

JJucxdxdy = (q1) j ui

ci+1, j - ci

+1/2,j

+ (q )i

ci, j- ci-1, j

(6)

где ы¡

+1/2, j +1, j

= (ui+1,j + ui,j )/2.

Вычислим интеграл, стоящий в правой части выражения (3):

Ц (Ц )'х dxdy = Ц (цсх )'х dxdy + Ц (цсх )'х dxdy.

В последнем выражении для определенности положим, что > , выделим из области Ц фрагмент Ц2, смежный с областью причем Ба = Ба (рис. 3).

Рис. 3. Схема заполненности областей

Л )XdxdУ = Л (цсх Ух^у + Л )Xdxdy =

^0 2 2

= ((41)г,. -(42 ),.)Я(ц)'хйМу + (к2).. Ц(цсх)\дхйу .

А А

Вычислим интеграл по области Б0:

у;+1/2 хг+1/2 у.+1/2

Я(Ц)'хйхйу = | йу | (цц)\йх = | цсх йу .

А0 у1-У2 хг-1/2 —1/2

Введем замену Ш = цсх и вычислим +1/2

У W

x +1

J — dx = J cXdx

и x

или Wi+1/2

(x+1 J Л 1

r dx /_ .

J — (+1 - ci)- U+ J и

\ x ^ У

x

и

В результате получим

Я (РР Ух^у-

Р

с+, . - с. .

'+1. j '' j

+1/2.у

Н

Р

с . -с. , .

' ' j '-1' j

-1/2'у

Н

Вычислим интеграл от функции (Рс )Х по области Д :

уу+1/2 Х'+1/2 уу+1/2

Я (РсХ )ХdxdУ = | ^ | (РсХ )Хdx = | рсХ 1Х'+"2 dУ:

Р

с'+1' .-с'' у

+1/2' у

Н

-р. (и +& )

Интеграл' стоящий в правой части выражения (3)' равен

(

' (ч2 )''у Р-1/2'у '

Л (РсХ — (Ч1)'' у Р' + 1/2' У 1 \ ± 2 /' у * ' 1/ 2'у ^

О) V НХ ' НХ

-((Ч1 )''у -(Ч2 )'' у )Р'' у (('у + АХ )) Ну .

В случае' если S0h > 5О ' результат будет аналогичным.

Подставим в уравнение (4) выражения (5)-(7) в результате чего получим:

(Чо)'.НхНу + (Ч1)'. и+1/2' у Ну +

(7)

+ (Ч2 )'

и.

—-Н + (ч3 ) V. .+,„ С-^Н +

о у /'' у ''у+1/2 0 Х

+ (Ч4 )''у^' (

с .-с . ,

'' у '' у-1

Н =

с'+1'у-с •■

(Ч1 )''у Р+1/2'у '+'Н ''у -(Ч2 )''у Р- 1/2'у ~ (Ч1 )''у -(Ч2 )''у | Р1'у ((у + РХ )) Ну +

с . - с. , .

' у ' -1' у

с. .,,-с.

(Ч3 )''у Р' у+1/2 "''У+1 "''У ( Ч4 )'' у Р''у-1/2

с .-с . ,

'' у '' у-1

(Ч3 )'' у -(Ч4 )'' у| Р'' у ('' у + Ау )) НХ +(Чо )'' у fi' у НХ Ну .

(8)

Разделим полученное выражение (8) на площадь ячейки НХНу' в результате

чего получим дискретный аналог уравнения диффузии-конвекции-реакции (1) с граничными условиями третьего рода (2):

с''у-с •■

(Чо)'' + (Ч1)''

с+, . - с. .

' + 1' у '' у

2Н

+ (Ч2)'

и..

с - с с - с

''У 1'У +(Ч3 )' у V''у.+1/2 +

2Н

2Н„

о.

+ («4 Ji, j—1/2

ci, j - cu-1

2h

= (q) i к

c+i,j - ci, j

+1/2, j

c. . — c. , .

■ j i—i, j

(q2 ),j K<-1/2,j ^2

(q), j —(q),

к, -

a c. . + ß

x i, J r'x

к

+ («з )i j K,

ci, j+1 — c

j+1/2

К

- — («4 )i j К,

ci, j—ci,j—1

j—1/2

к2

(«3 )i, j — («4 )i,j

aci,-+ß ( ч ,

к,- y к y +(q°)i,jyi,-.

(9)

Таким образом, получили дискретные аналоги операторов конвективного и диффузионного переноса в случае частичной заполненности ячеек.

Аппроксимация операторов диффузии-конвекции. Рассмотрим дискретные

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

е. -е. , .

' (10)

(«0 )i j uc'x = («1 )i j ui+1/2, j 1+1^;J , l,J +(«2 )i j ui —1/2, j — 1 —^J

2h

(q )i, - (кХ)'х - («1 )i, - К

+1/2, j

^j1- — ^ )i,j К i,J

2hx

c.. ,■ — c.

1,j

—1/2,j

h

(11)

(«1 )i, j — («2 )i,j

К

axci, j + ß h

где ц0, д1, ц2 - заполненности контрольных областей; а,Р - коэффициенты, стоящие в граничных условиях.

Для определения погрешности аппроксимации выражений (10), (11) будем доопределять расчетную область. Выражение (10) можно рассмотреть в случае (ч1. = (<?2. = 1, при этом можно утверждать, что погрешность аппроксимации

полученного выражения равна погрешности исходного выражения. Влияние граничных условий учитывается в выражении (11). Для определения погрешности аппроксимации выражения (11) нужно рассмотреть два случая: первый случай не учитывает влияния границы (ц1) = (д2). . = 1, второй случай учитывает

). . = 1; (д2 ) = 0 , так как аппроксимация (11) может быть записана через линейную комбинацию аппроксимаций, полученных в этих двух случаях. Таким образом, получим, что для определения погрешностей аппроксимаций выражений (10), (11) достаточно исследовать на точность следующие аппроксимации:

- дискретный аналог оператора конвективного переноса в случае отсутствия влияния границы области

ucr

u

ci+1, j ci, j

i+1/2, j

2h

u

ci, j ci—1, j

i —1/2, j

2h

(12)

- дискретный аналог оператора диффузионного переноса в случае отсутствия влияния границы области

ИХ - К

ci+1, j ci, j

i+1/2,j

-к

ci, j ci—1, j

i —1/2^.7

(13)

- дискретный аналог оператора диффузионного переноса в случае наличия влияния границы области

/ /\Х ,„ с1+1'у - с1,у аХс1,у + Рх

(РсХ ) х / 2 — Р+1/2'у -- Р'''У -Н-. (14)

Х

Погрешность аппроксимации оператора диффузии-конвекции. Найдем погрешность аппроксимации выражения (12)' для этого воспользуемся разложением в ряд Тейлора относительно узла (''' у) значения функций в узлах ('' +1' у) и

('• -1' у)

с'+1' у = с'' у +(с' у )х + ( у ) Н2 + 0 (Нх3 )' (15)

х х Н 2

с'-1'у = с''у - (с''у ) Нх + (с''у ) НТ + 0 (Нх3 ) . (16)

Подставим выражения (15)' (16) в аппроксимацию (12)' в результате чего получим:

/ и'+1/2 у + и'-1/2 У ' и'+1" ; - и'

ис — —

2 у ''" 4

Принимая во внимание выражения

' у/ V '+1/2' у '-1/2' у/ Г, У (с''У ) +-—А-- (с''У ) Нх + 0 ( ).

Щ+1/2'у + Щ-1/2'у = 2и'''у + 0 (НХ ) ' Щ+1/2'у - Щ-1/2'у = 0 (Нх )'

получим

и+,/9 . с'+''У с''У + и' ,„ . с''у с-и = и. . (с' . )' + 0(Н2). (17)

'+1/2'у 2Нх '-1/2'у 2Нх "'А ''у/ Vх/ v 7

Найдем погрешность аппроксимации выражения (13)' для этого воспользуемся разложением в ряд Тейлора относительно узла (''' у) значений функции с в

узлах ('' +1' у) и ('' -1' у):

с'+1'У = си +(си )Ч +(си ) Н2 + (с'' у) ^ + 0 (Нх4)' (18)

с'-1'У = с,у -(с''У )'Нх +(с''У ) ^-(с''У ) ^ + 0 (Нх4 ). (19)

Подставим выражения (18)' (19) в аппроксимацию (12)' в результате чего получим:

(РсХ )) — Р'+1/2'У-Р'-1/2'У (с'' У ) + Р'+1/2'У 2 Р'-1/2'У (с'' у ) + + ('+1/2'У -Р' -1/2'У )(' у + 0 (Нх2 ).

Принимая во внимание выражения

Р1+1/2'у + Р1 -1/2'у = 2Р'''у + 0 (НХ ) ' Р+1/2'у - Р-1/2'у = Ру ) Нх + 0 (% ) '

получим

Р-Р-1/2'У .Г-. = (р''У (с. ) | + 0 (Нх2). (2о)

Погрешность аппроксимации граничных условий третьего рода. Найдем погрешность аппроксимации выражения (14). Следует отметить' что слагаемое

Р у (ахс1 у + Рх) / Нх описывает влияние границ расчетной области и не вносит вклад в

'' У \ х '' у

погрешность аппроксимации [7]. Следовательно' его можно исключить при исследовании погрешности аппроксимации разностной схемы (14). Рассмотрим случай' когда ах = /Зх = о. В данном случае выражение (14) примет следующий вид:

с - с

( /-1 '+1' У '' У

Рх )х — 2Р' +1/2'у -. (21)

Тогда аппроксимацию уравнения диффузии в случае граничного узла можно записать в следующем виде:

с1' у - с' у = с+1' у - с1, у

2Н, ~Р'+1/2'у Н2 ' (22)

t х

где с - значение функции на предыдущем временном слое; с - значение функции на текущем временном слое; с = стс + (1 -о)с .

Воспользуемся разложением в ряд Тейлора относительно узла (''' у) значения функции в соседнем узле

Н:

Е/ \{m) nx

(J) . (23)

Tl. I z_t \ i. II .

m!

Введем фиктивный узел ('' -1' У). Разложение функции с в ряд Тейлора в данном узле относительно узла ('\ у) запишется в следующем виде:

Е, \(т! / ,\т Нх

(с'''.) (-1) т^. (24)

т=о т:

Можно доопределить задачу в нерасчетной области следующим уравнением:

сц - сц сц - с -1 у

—-У = -Р'-1/2 у 'у 2 'у . (25)

2Нt 1/2'у Н2

t х

Нетрудно видеть' что если сложить левые и правые части выражений (22)' (25)' получим разностную схему' аппроксимирующую уравнение диффузии со вторым порядком точности. Из выражений (22)' (25) в случае симметричности коэффициента

диффузии относительно границы расчетной области (Р -1/2' у = Р+1/2' у) следует равенство с'+1' у = с'-1' у. Приравняв правые части выражений (23)' (24)' получим:

Р^^. (26)

Подставим (23) в (22) с учетом полученного выражения (26)' получим:

c — c

^.J . J

h

■■MtTi/2.J ((.J)' + О(hx2).

В силу симметричности коэффициента диффузии относительно границы расчетной области имеет место выражение

(Р'у (у ^ = 2Р'+1/2'у с'+1'^ с''У + 0(Нх2).

Таким образом, получили, что погрешность аппроксимации разностных схем, аппроксимирующих операторы диффузионного и конвективного переноса (10), (11) равна 0(к2х ).

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

Ь( с( Р )) = А( Р )с( Р)- X В(р, 6 )с( 6 ) = р(р), (27)

2еШ' (Р)

где Ь - некоторый сеточный оператор; Р = ( л1, yj) - центр шаблона;

Ш '( Р )= {б! ( X; + 1 , Уj ) , й2 ( X;-! , Уj ) , й ( X; , Уj ц ) , 64 ( X; , Уj^ )} - окрестностЬ центра

шаблона.

Коэффициенты сеточных уравнений, стоящие в окрестности центра шаблона, примут вид

В( Р, й ) = а В( Р, 62 ) = ^ В( Р, 6з ) = ^ В( Р, 64 ) = ^

и(61) + и(Р) | м(61 ) + ^( Р) 4к 2к2

х х

и( 62) + и( Р) + ц(С62 ) + ^( Р)

4к

2к2

(

К 63) + К Р) + м( 63 ) + ^( Р)

4кУ 2к2

' у( 64) + у( Р) + м( 64 ) + ^( Р)А

4к

2к~

91 ( Р),

92 ( Р ), 9з( Р ) 94 ( Р ).

У У

Коэффициент, стоящий в центре шаблона, запишется так:

А( Р ) = + В( Р, 6,) + В( Р, 62) + В( Р, 63) + В( Р, 64) +

+а

Ы Р)-92 (Р )| кг + 9 ( Р)-94 ( Р )| ^ Ц(Р) .

у /

Для записи правой части сеточного уравнения нам понадобятся вспомогательные коэффициенты

В( Р, 66 ) = ( 1 В( Р, 67 ) = (1

и( 61) + и( Р) + 61 ) + м( Р)

V 4кх + 2к2 ,

91 (Р),

и( 62) + и( Р) + 62 ) + м( Р)

4к

2к

92 (Р )

В( Р, 68 ) = ( 1

К 6з) + у( Р ) + ^( 6з ) + м( Р)

4к„

2к2

9з ( Р )

В( Р, 69 ) = ( 1

у( 64) + у( Р ) + ^( 64 ) + М Р)

4к

2к2

94 (Р)

B (P, 05 )= - B (P, Q6)-B (P, Q7)-B (P, ß8)-B (P, 09) +

(

,Г

л

+ (1 (Р)-^(Р))-^ + \чз (Р)-44(Р^т2- р(Р)-

I ^ ку ,1

При этом правые части сеточных уравнений примут вид Р(Р) = 4о (Р)/(Р)-к (Р)-92 (Р)|р(Р-

-|4з (Р)-44 (Р)ИР+ В(Р,Й5)(Р) + В(Р,&)) +

+в (р, е7) (е2)+в (р, е8 )с &)+в (р, & )с (&).

Достаточные условия устойчивости. Условием применимости принципа максимума являются:

Л (Р)> 0, В (Р, б )> 0, О (Р) = Л (Р)- £ В (Р, б)> 0. (28)

ОеШ '(Р)

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

B (P, Q, ) = а

u (Q) + u (P) ^(Q ) + ^(P)

4h

2h2

q, (P)> 0, i = 1,4 :

следует

hx < 2minГм/|u|) < 2 Qi) + P). x V 1 V I u(Qi) + u(P)

Ограничения на шаг по времени имеют вид

' (1 >)) + )) + (1 • т< -(1 -а)--+(1 -а

2min((P))

(29)

(1 -а)

qo (P

г/ Г

|qi (p) -q2 (P)\r + |q3 (P)-q4 (P)!^ M(P)

к

y

В случае неявной схемы (а = 1) отсутствует ограничение на шаг по временной координате. Условием монотонности явной схемы для уравнения диффузии в случае граничных условий первого и второго рода является: т < tf / (2min(m(P))).

Проверка устойчивости. Согласно принципу максимума выполнения условий (3), имеет место следующая оценка:

IHI c lF / DL, (30)

где F - значение функции, стоящей в правой части сеточного уравнения; D -диагональное преобладание.

Следовательно, для исходного уравнения имеет место оценка:

< с0 +т

i\fk

qi (p )-q2 (P )\y+q (P )-q4 (P

|qi (P)-q2 (P+ q (P)-q4 (P)\Г

Из полученного выражения следует устойчивость исследуемых математических моделей волновых процессов по начальным данным, граничным условиям и правой части.

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

Проверим баланс вещества. Для этого просуммируем уравнение (4) по всей расчетной области с учетом е1+1/2, ■ = ( е1+1, у + е1, ■) / 2 , в результате чего получим

X ( \jCij = X ( 4о\jCij + X т( 4о );,у( их + );Ау +

;е[1, -1], ;Е[1, -1], ;Е[1, -1],

j■E[l,Му -1] j■E[l,Му -1] j■E[l,Му -1]

e[1, Nx-jefl. Ny-1|

a*Ci,i + ß.l/ 4 / Ч I.. ayCi,j + ßy

- hX(i r | ( qi \j -();j| ßi^h—1+|( L -();j| ß h

•e ' " y у

+

+г X ( qo)jfi,j. (3i)

;e[i, nx -il,

jE[l,Ny -l]

где ux. vy - операторы первой разностной производной.

Из полученного равенства следует. что оператор диффузии не вносит вклад в суммарную концентрацию вещества. Оператор конвективного переноса не является консервативным в случае невыполнения условия несжимаемости среды

ux + vy = 0 . Если среда сжимается (расширяется). то концентрация увеличивается

(уменьшается). что иногда удобно при моделировании реальных физических процессов.

Проинтегрируем данное уравнение по расчетной области G :

И ( < + ucX + vcy )dS = Ц i (ficX )'x + ( Mc'y)y + f ) dS . (32)

G G V у

Применим теорему Гаусса - Остроградского

JJc'dS = JJ ( c (ux + vy) + f )dS + ф ( -Vc + ßgradc. n0dl). (33)

G G у

Выражение (26) с учетом (27) может быть записано в следующем виде:

JJC'dS = JJ( c( ux + Vy) + f)dS + ф( -(un0 + vn0)c + ß(anc + ß))l . (34)

G G у

Нетрудно видеть. что уравнение (31) является дискретным аналогом уравнения (34).

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК 1. Самарский А.А. Теория разностных схем. - М.: Наука, 1989.

2. Сухинов А.И., Чистяков А.Е., Бондаренко Ю.С. Оценка погрешности решения уравнения диффузии на основе схем с весами // Известия ЮФУ. Технические науки.- 2011. - № 8 (121). - С. 6-13.

3. Сухинов А.И., Тимофеева Е.Ф. Чистяков А.Е. Построение и исследование дискретной математической модели расчета прибрежных волновых процессов // Известия ЮФУ. Технические науки. - 2011. - № 8 (121). - С. 22-32.

4. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978.

5. Коновалов А.Н. К теории попеременно-треугольного итерационного метода // Сибирский математический журнал. - 2002. - № 43:3. - С. 552-572.

6. Чистяков А.Е. Теоретические оценки ускорения и эффективности параллельной реализации ПТМ скорейшего спуска // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 237-249.

7. Чистяков А.Е. Об аппроксимации граничных условий трехмерной модели движения водной среды // Известия ЮФУ. Технические науки. - 2010. - № 6 (107). - С. 66-77.

8. Фоменко Н.А. Моделирование гидродинамических процессов при обтекании корпуса судна // Известия ЮФУ. Технические науки. - 2011. - № 8 (121). - С. 139-147.

9. Сухинов А.И., Чистяков А.Е., Проценко Е.А. Построение дискретной двумерной математической модели транспорта наносов // Известия ЮФУ. Технические науки. - 2011. - № 8 (121). - С. 32-44.

Статью рекомендовал к опубликованию д.т.н., профессор Я.Е. Ромм.

Сухинов Александр Иванович - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: sukhinov@gmail.com; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 89281021106; д.ф.-м.н.; профессор.

Чистяков Александр Евгеньевич - e-mail: cheese_05@mail.ru; тел.: 88634371606; кафедра высшей математики; доцент.

Фоменко Наталья Алексеевна - e-mail: fomenko.n86@mail.ru; тел.: 89034855580; кафедра высшей математики; ассистент.

Sukhinov Alexander Ivanovich - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: sukhinov@gmail.com; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +79281021106; dr. of phis.-math. sc.; professor.

Chistyakov Alexander Evgenjevich - e-mail: cheese_05@mail.ru; phone: +78634371606; the department of higher mathematics; associate professor.

Fomenko Natalya Alexeevna - e-mail: fomenko.n86@mail.ru; phone: +79034855580; the department of higher mathematics; assistant.

УДК 532.5.031

А.И. Сухинов, А.Е. Чистяков, М.Д. Чекина

ОПИСАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ПРОЦЕССА ПЕРЕМЕЩЕНИЯ СЫПУЧИХ ВЕЩЕСТВ С ИСПОЛЬЗОВАНИЕМ УРАВНЕНИЯ СЕН-ВЕНАНА

Рассматривается описание математической модели процесса перемещения сыпучих веществ с использованием уравнения Сен-Венана. Для системы дифференциальных уравнений, описывающих данную модель, получены дискретные аналоги. С помощью метода контрольных объемов, применяющегося для дискретизации модели, достигнута высокая точность аппроксимации, в том числе и на грубых сетках. На основе полученной модели разработан комплекс программ, с помощью которого производится численное моделирование