Численное моделирование плоских течений неньютоновской жидкости со свободной поверхностью непрямым методом граничных элементов Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 11, № 5, 2006

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПЛОСКИХ ТЕЧЕНИЙ НЕНЬЮТОНОВСКОЙ ЖИДКОСТИ СО СВОБОДНОЙ ПОВЕРХНОСТЬЮ НЕПРЯМЫМ МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ*

В.А. ЯКУТЕНОК, М.Н. ШТОКОЛОВА Томский государственный университет, Россия e-mail: [email protected]

A numerical algorithm for modeling of nonlinear viscous liquid flows is presented. It is based on application of the boundary elements method. The algorithm was tested on the problem of flow in the plane channel, the computation accuracy has been verified by comparing the present numerical results with the analytical solution.

Введение

Задача о вращающемся течении внутри полого горизонтального цилиндра интересует исследователей на протяжении многих лет вследствие широкого спектра ее приложений для нужд промышленности. Математическое моделирование процесса перемешивания полимерных композиций в объемном смесителе сопряжено со значительными трудностями. Эти трудности связаны, прежде всего, с пространственным характером процесса перемешивания, неньютоновским характером поведения перемешиваемой среды и наличием свободной поверхности. Для понимания механизма перемешивания, качественного представления о степени влияния определяющих факторов используют различные упрощения математической формулировки задачи. Основным упрощением является предположение о том, что смеситель представляет собой бесконечно длинный горизонтальный цилиндр, частично заполненный высоковязкой жидкостью и вращающийся вокруг собственной оси. Для такого случая в предположении ньютоновского характера течения известны работы как зарубежных [1, 2], так и отечественных авторов [3]. В то же время, указанные исследования не затрагивают вопросов о влиянии отклонения реологии среды от закона Ньютона. Известны несколько попыток аналитического исследования пленочных течений неньютоновской жидкости внутри горизонтального цилиндра, где за основу взяты степенная реологическая модель [4, 5], модель Бингама [6], а также обобщенная модель ньютоновской жидкости [7]. Однако в упомянутых работах рассматривается лишь пленочный режим течения, а также не представляется возможным количественное определение кинематических и динамических характеристик внутри перемешиваемой жидкости. В связи с этим

* Работа выполнена при финансовой поддержке Министерства образования и науки РФ и CRDF в рамках программы BRHE (грант № 016-02) и Российского фонда фундаментальных исследований (грант № 06-08-00107-а).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

Сделанное предположение о том, что течение в смесителе качественно совпадает с течением неньютоновской жидкости в горизонтальном вращающемся частично заполненном бесконечном цилиндре, приводит к необходимости решения задачи из класса нелинейных плоских задач гидродинамики ползущих течений со свободной поверхностью. Допущение о применимости ползущего приближения основано на малости характерных чисел Рейнольдса течений, реализующихся в смесителе. Для течений, реализуемых в объемных смесителях, максимальное значение И,е = 0.71. С другой стороны, хорошо известно, что условие применимости ползущего приближения И,е < 2. Для решения нелинейных плоских задач гидродинамики течений со свободной поверхностью можно использовать различные конечно-разностные методы [8]. При рассмотрении течений в областях сложной геометрии, с большими деформациями свободной границы, при наличии нескольких, значительно различающихся характерных размеров возникают вычислительные трудности, связанные с построением конечно-разностных сеток, которые, кроме того, должны адаптироваться к изменениям свободной поверхности. В связи с этим актуальность приобретают поиски подходов, позволяющих упростить и унифицировать алгоритмическую процедуру численного решения. Как известно, одним из путей преодоления этих трудностей является переход к системе граничных интегральных уравнений, эквивалентной исходной системе дифференциальных уравнений [9, 10]. Учет реологического поведения, отличного от модели Ньютона, приводит к дополнительным сложностям, способам преодоления которых посвящена данная статья.

1. Математическая постановка задачи о течении степенной жидкости

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

Тц = 2кАп ^ё], (1)

где Т] — компоненты девиатора тензора напряжений; к — показатель консистенции; п — индекс течения (показатель нелинейности); ё] = (дгг/дх] + дг]/дхг)/2 — компоненты тензора скоростей деформаций; гг — компоненты вектора скорости; хг — декартовы координаты; А — второй инвариант тензора скоростей деформаций (интенсивность скоростей деформаций) [11]. Уравнение (1) (так называемая модель Оствальда де-Вилля) описывает поведение псевдопластичных (п < 1) и дилатантных (п > 1) жидкостей. При п =1 имеем ньютоновскую жидкость, при этом коэффициент консистенции к совпадает с коэффициентом динамической вязкости. При наличии предела текучести модель Оствальда де-Вилля обобщается реологическим уравнением Шульмана [11]. С методологической точки зрения такое обобщение не вносит принципиальных изменений в построение математической модели и метода решения.

В приближении ползущего течения в присутствии силы тяжести течение жидкости описывается уравнениями Стокса, которые можно записать в виде

да7 , ,

3X7 + Р9г = 0, г,3 = 1 2, (2)

где а7 = — рбу + ту — компоненты тензора напряжений; р — давление; р — плотность; д. — компоненты вектора ускорения силы тяжести.

К уравнению (2) следует присоединить уравнение неразрывности

ду . .

зх = °- (3)

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

у = V(x), X € Б\, ¿¿(х) = а у Пу = 0, X € Б2, (4)

где Б2 — свободная поверхность; Б\ = Б П Б2; пу — компоненты внешней нормали к Б2.

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

дf

^ + г&гаЛ/ = 0, (5)

где /(х, ¿) — функция, задающая Б2. В лагранжевой форме кинематическое условие имеет вид

у=!■ (б)

Условия (5), (6) используются для нахождения формы свободной границы. Выбор варианта определяется характером поведения свободной поверхности для наиболее точного вычисления ее формы. Реологическое уравнение (1) представим в виде

ту = 2Веу, (7)

где В = кЛп-1 — коэффициент эффективной вязкости. Для перехода к безразмерным переменным введем следующие масштабы: для координат — Я (характерный размер области, занятой жидкостью), для скорости — V, для напряжений — к(У/Я)п. Тогда (сохраняя прежние обозначения для безразмерных переменных) уравнение (2) примет вид

3а7 + - = 0, г,з = 1, 2. (8)

Здесь Uij = -p8ij + Tij — безразмерные компоненты тензора напряжений; Tij = 2Bej — безразмерные компоненты девиатора тензора напряжений; B = An-1 — безразмерный коэффициент эффективной вязкости; Wi = W(gei)/\g\, W = (p\g\R/k) (R/V)n, ei — орты декартовой системы координат.

Следовательно, кроме геометрических, имеется всего один определяющий параметр — число W, которое, являясь отношением чисел Рейнольдса и Фруда, модифицированных для нелинейно-вязкой жидкости, характеризует соотношение гравитационных и вязких сил в потоке жидкости. Для последующей линеаризации преобразуем уравнение (8) следующим тождественным образом:

о о

dj + Wi = dj [2(1 - B)eij ] ' ij = !' 2' (9)

где aij = -pSij + 2eij — линейная часть тензора напряжений. В случае B = 1 будем иметь линейное уравнение, описывающее течение ньютоновской жидкости. Следовательно, в таком представлении правая часть уравнения (9) описывает нелинейное поведение жидкости, обусловленное зависимостью вязкости от скорости сдвига.

В уравнении (9) можно избавиться от постоянного слагаемого Wi, вводя потенциал для силы тяжести: ф = WiXj. Тогда, сохраняя прежние обозначения, имеем

о о

dj = dj12(1 - B)eij]' = 1'2' (10)

но теперь aij = — (p — ф)5^ + 2ej. Далее под величиной давления будем понимать его модифицированное значение p — ф, однако для удобства обозначение для модифицированного aij оставим прежним. Соответственно видоизменится граничное условие на свободной поверхности, которое следует использовать в виде

ti(x) = WjXjni, x £ S2, (11)

где ti(x) = aijnj, а под aij понимаются безразмерные компоненты полного тензора напряжений. Таким образом, математическая формулировка задачи включает уравнения (3), (10) с граничными условиями (4)-(6) с учетом (11). В качестве начального условия необходимо задавать форму свободной границы в начальный момент времени.

2. Граничные интегральные уравнения для нелинейной задачи

Для перехода к гранично-интегральной формулировке задачи запишем уравнение (10) в виде

§х?=*'• (12)

д

где Фг = —— [2(1-В)е]] — нелинейная векторная функция, которую будем рассматривать дх]

как плотность источников, распределенных по области течения.

Как известно, линеаризованная система уравнений Стокса для случая единичной сосредоточенной силы, действующей вдоль координатной оси ] и приложенной в точке £, имеет вид

АУ, — grad р, = 5(х — £)е,, divVj = 0, (13)

где 5 — дельта-функция Дирака. Дифференцирование проводится по переменной х, £ рассматривается как параметр. Сингулярное решение системы (12) записывается следующим образом [12]:

^ = — 4п

5 п 1 (хг — £г)(х, — £,) °гз |х — £| + |х — £ I2

£ I

= 1 1

р = 1п £|'

(14)

Далее будем использовать обозначения: г = |х — £|, уг = хг — £г, 5г, — символ Кронекера. Подставляя (14) в выражение для напряжений, можно получить

Уг Ук у,

°гк3 = ' (15)

Следовательно, для усилий в точке х поверхности с нормалью Пк(х) имеем

УгУ, Ук Пк , .

тг, = яг, Пк = —ПГ4—' (16)

Соотношения (14), (16) позволяют перейти к интегральной формулировке поставленной задачи. Для этого реальную область течения заменим фиктивной, отличающейся тем, что она считается помещенной в неограниченную двумерную область, заполненную ньютоновской жидкостью. На границе фиктивной области будем считать распределенными сосредоточенные силы неизвестной интенсивности ф, (£). Внутри области течения дополнительно распределены внутренние источники, интенсивность которых определяется функцией Ф,(г). Здесь введено обозначение г для переменной интегрирования по области течения П, аналогичное обозначению £ для переменной интегрирования по границе Б. Тогда скорости уг и усилия ¿г в любой внутренней точке х можно найти сверткой фундаментальных сингулярных решений (14), (16) с функциями ф, (£), Ф, (г) по следующим формулам:

Уг(х) = I Уг,(х,£)ф,(£)^Б(£) + [ Уг,(х,г)Ф,(г)^(г); (17)

*г(х) = у из(х,£)ф,(£)^Б(£) + у из(х, *)Ф,(г^ВД. (18)

5 П

Уравнения (17), (18) удовлетворяют дифференциальному уравнению (12) в любой внутренней точке и определяют решение при выполнении граничных условий. Функция плотности распределения внутренних источников Ф, (*) заранее неизвестна и может быть вычислена только для известного поля скорости. Это приводит к необходимости аппроксимации Ф,(*) (например, с использованием поля скорости для вязкой жидкости) с целью дальнейшего уточнения решения при помощи того или иного итерационного процесса. Особенности во вторых интегралах уравнений (17), (18) при * = х являются интегрируемыми по области решения.

1

Последний шаг в получении интегральных граничных уравнений состоит в помещении точек наблюдения на границу Я (х ^ х0, х0 Е Я). В этом случае интегральные уравнения (17), (18) становятся сингулярными, так как фундаментальные решения не ограничены в точках х = £. Первый интеграл, имея в ядре логарифмическую особенность, может быть вычислен в обычном смысле. Второй интеграл в подынтегральном выражении имеет особенность 1/г, поэтому его следует понимать в смысле главного значения по Коши. Свободный член вычисляется стандартным методом теории потенциала, и в точке хо, в которой существует единственная касательная к границе, можно записать

¿¿(х) = 1 фг(хо) + 1 из(х,£)ф,(£)^(£) + / ¿ц(х, г)Ф, (г)ЛВД. (19)

5 П

Граничные интегральные уравнения (17) и (19) при х ^ х0 и заданных в соответствии с (4), (11) уг(х0) и ¿¿(хо) позволяют получить значения неизвестных граничных сил ф(£). Это возможно сделать, только если известна функция Фц(г). По этой причине возникает необходимость линеаризации уравнений (17), (19) и организации итерационного процесса.

3. Численный метод решения

Найти аналитическое решение нелинейной системы интегральных уравнений (17), (19) не представляется возможным даже для областей простой формы, поэтому она решается численно. Разобьем границу области, занятой жидкостью, на N прямолинейных отрезков (элементов) и будем считать функцию фj (£) постоянной в пределах элемента. Область П, занятую жидкостью, разделим на М выпуклых четырехугольников, внутри каждого из которых функцию Фj(г) примем постоянной. Тогда уравнения (17), (19) в дискретизиро-ванной форме приобретут вид

N М

^(х?) = £фj(£*) ] Уц(хР,е9)^(£9) + $>,(г1) ] Уц(х?^)^); (20)

9=1 АЯч 1=1 АП1

N М

¿¿(*0) = £фц(£9) ] ¿ц(х0,е9)^(£9) + $>,(г1) ] ¿ц(хО,^)^1). (21)

9=1 АЯч 1=1 АПг

Здесь хр — некоторая точка элемента с номером р, которую будем называть узлом. В качестве узлов будем выбирать середины граничных элементов. Из системы линейных алгебраических уравнений (20), (21) берутся N уравнений, причем уравнения (21) соответствуют элементам на свободной поверхности; на всех других элементах заданы скорости, поэтому для них используются уравнения (20).

Интегралы У Уц(х0,£9(£9) и J ¿ц(х0,£9(£9) в уравнениях (20), (21), являющи-

АЯч АЯч

еся коэффициентами системы линейных уравнений, можно вычислить аналитически. Для

этого вводится декартова система координат уг (г = 1, 2), нормально связанная с элементом д (не путать с ранее введенным обозначением уг = хг — £г). Причем ось у1 направлена по внешней нормали к элементу, ось у2 — в положительном направлении касательной (против часовой стрелки). Начало системы координат уг находится в точке р, т. е. уг = ец(хц — хр).

Здесь в, — направляющие косинусы осей. В этом случае, если к тому же использовать полярные координаты г, в, (рис. 1) можно получить

..... , _

4п

Д5 в

V??¿Б(£?) = —<1 (1п г - 1)У2 - У1 (в,-в2г + в2,в14) 1п Г - в2,в2г (у2 - 2^16) ^ ,

' ?

¿Б(£?) = в,ви [те 1 (6 + Бт 6 еоэ 6) + п2Бт2 6 +

2п

Д5 в

(22)

+ (в1, в2^ + в2,в1^) [п 1 й1п2 6 + п2 (6 - бш 6 еоБ 6)] +

+ в2, в2^ [п 1 (6 - б1п6соб6) - п2 (бш2 6 - 21п г)]

?+1

?

где щ = в, п,.

Уравнения (22) верны для всех узлов х?, за исключением случая д = р, т. е. совпадения точек наблюдения и приложения нагрузки. В этом случае

6? ^ -п/2, 6?+1 ^ п/2, П1 = 1, П2 = 0,

г? = г?+1 ^ ДБ?/2, у? = у?+1 ^ 0.

Тогда запишем

¿Б(^) = Д^ [1п(ДБ?/2) - 1] - в2,в2г ¡>,

А£ в

1 (23)

) %^(С) = ^,

А£ в

где ДБ? — длина элемента д.

Второе уравнение (23) показывает, что вклад в интеграл при д = р дает лишь свободный член, главное же значение по Коши на элементе равно нулю. Используя формулы (22), (23), можно получить систему уравнений, которая в матричной форме будет иметь вид

2Ж х1

' / V?? ¿Б / ¿Б

1_А£

2Ж х2Ж

^ + • (24)

2Ж х 1 2Ж х 1

Система (24) решается стандартным методом Гаусса с выбором главного элемента, и в результате определяется вектор-столбец {ф9}. Для учета нелинейного слагаемого {ФР} следует организовать итерационный процесс с использованием, например, метода простой

итерации. Интегралы У Уц(хР,гг)^П(гг), У ¿ц(хР,гг)^П(гг) можно вычислять с помощью

АП1 АП1

численных квадратур Гаусса. Скорости, давление и напряжения рассчитываются в любой точке области с помощью численного интегрирования по дискретным аналогам соотношений, подобных уравнениям (20), (21).

После нахождения значений скоростей в узлах, принадлежащих свободной границе, рассчитывается ее новое положение с помощью кинематического условия, записанного в виде (5) или (6). Для реализации условия (5) можно использовать разностную схему с разностями против потока, которая в случае использования декартовых координат имеет вид

/п+1 _ /п Д / п

Лг Лг _ п п Д/г

Уо. — У

Дг 24 14 /пДх1' Д/Г = {/Г_ /Г-1,УП14 > 0, (25)

/ _ /Г,уПг < 0.

Здесь п — номер шага по времени; Дх1 — шаг разностной сетки по координате х1; г — номер узла на свободной поверхности; Дг — шаг по времени. Условие устойчивости схемы (25) имеет вид

Дхг

At < min

г

Граничное условие (6) дает следующие разностные формулы для вычисления новой формы свободной границы (схема Эйлера):

= + Atvn, i = 1, 2. (26)

Вместо (26) возможно использование схем повышенного порядка точности (типа предиктор-корректор). В соответствии с (25) или (26) находятся новые координаты узлов, принадлежащих свободной границе. Далее организуется новый итерационный процесс для учета нелинейного слагаемого.

Таким образом, данную задачу можно характеризовать как квазистационарную, поскольку для заданной формы границы находится решение, используемое для нахождения новой формы. Если рассматривается стационарный случай, то в итоге должна получаться стационарная форма свободной поверхности. На первой итерации первого шага по времени должна решаться линейная задача, что позволит найти первое приближение для (Фр|. Далее в качестве первого приближения можно использовать значения с предыдущего шага по времени.

4. Решение тестовой задачи о течении нелинейно-вязкой жидкости в плоском канале

Выше изложен непрямой метод граничных элементов для решения задач о течении неньютоновской жидкости со свободной поверхностью. В таком варианте этот метод ранее не

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

эффективно вычислять интегралы / V,,(4,г')<Ш(;') и / )««(*-) по элементам

ДП1 ДП1

области решения. В работе [9] предлагается использование полуаналитической схемы интегрирования для элементов области ДП1, если € ДП1. В то же время представляется возможным применение квадратур Гаусса повышенного порядка точности с использованием внутренних точек области, так как, несмотря на особенности в подынтегральной функции, сами интегралы существуют в обычном смысле. При этом особенность должна находиться на границе ДП1. В этом случае алгоритм решения значительно упрощается ввиду единообразной процедуры вычисления рассматриваемых интегралов. Для этого необходимо решение тестовой задачи с известным аналитическим решением, что позволит также подтвердить достоверность результатов применения разработанного метода к численному моделированию течения неньютоновской жидкости в других задачах.

В качестве тестовой будем рассматривать задачу об установившемся течении нелинейно-вязкой жидкости в плоском канале. В качестве характерного размера выберем ширину канала Ь, в качестве характерной скорости — среднюю скорость течения V. Тогда можно получить следующее решение, описывающее поле течения [11]:

VI = 0, = 2П+1(1 - |2ж1 - 1|^), (27)

п +1

где VI — поперечная, а V2 — продольная составляющие вектора скорости; ж1 — поперечная координата с началом на оси канала.

Сформулируем соответствующую краевую задачу, точное решение которой выражается формулами (27). Область решения представляет собой прямоугольник с шириной I и высотой Н (в расчетах бралось Н = 1). Пусть на нижней (входной) границе в качестве краевого условия заданы компоненты скорости в соответствии с (27). На левой и правой границах (твердых стенках) задаются условия прилипания

V = 0, ^ = 0. (28) На верхней (выходной) границе зададим краевое условие в виде

V = 0, ¿2 = 0. (29)

Здесь ¿2 = °"22 = —Р — усилие в направлении оси ж2, действующее на выходной границе, равное давлению с обратным знаком. Без ограничения общности можно считать, что р = 0, так как давление — величина аддитивная.

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

Для вычисления интегралов У V,,(ж0,г)^П(г), У (жР,г)^П(г) разобьем область реп п

шения на прямоугольники ДП1 сеткой прямых, параллельных осям координат. В каждом прямоугольнике интегралы У V,,)^П(г), У (ж0,г)^П(г) будем вычислять по

ДП1 ДП1

36-точечным формулам Гаусса вида

Vj (xP, z)dQ(z)

дпг

íjj(xP, )

hih2

6 6

'y ] ^ ] A2 Amvíj (x0) z2 )

2=1 m=1

66

(30)

^ ^ ^ ^ A2 Ат^г7 (x0, zl , z2 )

ДПг

2=1 m=1

где Ak, Am — весовые коэффициенты; zl,z2 — точки Гаусса; h1;h2 — поперечный и продольный размеры прямоугольника ДП1 [9]. Формулы (30) будем единообразно использовать как для областей ДП1, не содержащих точку xp, так и для приграничных ячеек, у которых на одной из сторон находится эта точка.

В данном случае распределение плотности внутренних источников

Ф,;

д

dxj

[2(1 - B)eij]

(31)

можно вычислить аналитически с использованием (27). Тогда будем иметь

Ф1 = 0, Ф2

4(2n + 1)

n

2(2n + 1)

n

n— 1

|1 - 2x11

1 -n n

n

(32)

Выражения (32) вычисляются в центрах прямоугольников ДП1, умножаются на (30) и суммируются по всем /. В результате формируется система линейных алгебраических уравнений (24), решение которой позволяет по формулам (20) вычислить не известные заранее значения скоростей на выходной границе.

Для решения задачи составлена программа расчета. Результаты вычислений приведены в табл. 1. Здесь

max IV2 — V2q I

1<q<nb

E2

\

nb _

E(V? — V? )2

q=1_

nb

E3

max

1<q<nb

V? — V?

%

— отклонения приближенного решения от точного в нормах С и (максимальное и сред-неквадратическое), а также относительное отклонение; — точные значения скоростей в соответствии с (27); — вычисленные значения скоростей; пЬ — количество граничных элементов на выходной границе; пт = 4пЬ — общее количество элементов. Причем на всех четырех сторонах области течения количество элементов одинаково и равно пЬ. Коэффициент нелинейности п = 0.5. Результаты вычислений показывают, что разработанный метод обладает аппроксимационной сходимостью и позволяет получать решения с точностью менее 1 % уже при пт = 32.

Таблица 1. Результаты расчетов для коэффициента нелинейности n = 0.5

nm 8 16 32 64 128 256

E1 4.310-7 1.210-2 1.310-2 5.910-3 2.610—3 1.310-3

E2 4.310-7 1.110-2 8.210-3 3.610-3 1.510—3 6.010-4

E3, % 3.710-5 1.5 9.710-1 8.810-1 7.610—1 7.310-1

Таким образом, использование квадратур Гаусса для вычисления интегралов / „«.,)«(, ) и / „в задачах о течении ненью— жидкости

ДП1 ДП1

со свободной поверхностью вполне обоснованно.

Однако при решении реальных задач о течении нелинейно-вязких жидкостей воспользоваться заранее известным аналитическим выражением для распределения плотности внутренних источников, как в данном случае, не представляется возможным. Поэтому описанный выше численный алгоритм был также подвергнут проверке на аппроксимаци-онную сходимость при условии, что плотность внутренних источников Ф, первоначально неизвестна и определяется в ходе расчетов численно с использованием выражения (31). В таком случае организуется полный итерационный процесс в соответствии с методом простой итерации, как было отмечено выше. Плотность внутренних источников Ф, на первой итерации принимается равной нулю, и решается тестовая задача о течении в канале для ньютоновской жидкости. Далее по вычисленному полю скоростей и (31) вычисляются Ф, во всех центрах прямоугольников с использованием конечно-разностных разложений производных, входящих в (31), эти значения умножаются на выражения (30), и первоначальная система линейных алгебраических уравнений дополняется учетом внутренних нелинейных источников в правой части (24). Решение полученной таким образом системы дает новое распределение скоростей и, в свою очередь, новые Ф,. Процесс повторяется до установления Ф, и скоростей. Сходимость процесса определялась в ходе вычислений с точностью до е = 10-8. Решения, полученные в ходе численного эксперимента, сопоставлялись с аналитическими решениями для профиля скорости в диапазоне изменения коэффициента нелинейности п от 0.1 до 1.5. Сравнение с аналитическим решением выполнялось, как и прежде, в нормах Е1, Е2, Наилучшая аппроксимационная сходимость в нормах Е1,Е2 (Е1 ~ 10-3,Е2 ~ 10-4) достигалась при значениях коэффициента нелинейности в диапазоне 0.5 < п < 1.3, а в норме Е3 (Е3 ~ 1%) — при 0.8 < п < 1.3; при 0.5 < п < 0.7 Е3 = 1-5%. Однако при п < 0.5, а также при п > 1.3 итерационный процесс расходится. В расчетах количество итераций не превышало 50, в большинстве расчетов процесс сходимости завершался гораздо раньше.

В табл. 2 представлены результаты расчета для коэффициента нелинейности п = 1.2. Как и в предыдущем тесте, здесь пт = 4пЬ. Характер сходимости итерационного процесса иллюстрируется рис. 2-4.

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

Результаты решения тестовой задачи доказывают эффективность разработанного метода в широком диапазоне изменения коэффициента нелинейности 0.5 < п < 1.3. В то же

Таблица 2. Аппроксимационная сходимость итерационного процесса, п = 1.2

пЬ 5 9 17 33 65

Е1 2.0810-2 7.2910-3 2.6310-3 2.3110-3 1.6110-3

Е2 1.3610-2 3.9610-3 1.3710-3 6.5810-4 3.1810-4

Е3, % 1.42 0.48 0.83 1.38 1.90

1.6-1

К -

1.20.80.4-

0Н—'—|—'—|—'—|—'—|—'—| 0 0.2 0.4 0.6 0.8 х 1

Рис. 2. Профиль скорости У2(х\), п = 0.5, треугольниками отмечено аналитическое решение. 0.0250.02-

0.0150.010.005-

0

Рис. 3. Итерационная сходимость в норме Рис. 4. Профиль скорости V2(x\), n = 1.2, тре-E2, n =1.2, M — число итераций. угольниками отмечено аналитическое решение.

время необходимо отметить, что поведение полимерных масс, используемых в химической технологии, удовлетворительно описывается уравнением (1) с коэффициентом нелинейности 0.7 < n < 0.9. Это позволяет сделать вывод о возможности применения непрямого метода граничных элементов в представленном варианте для решения более сложных задач, в частности для моделирования течения неньютоновской среды в объемном смесителе.

Список литературы

[1] Orr F.M., Scriven L.E. Rimming flow: numerical simulation of steady viscous free surface flow with surface tension //J. Fluid Mech. 1978. Vol. 84. P. 145-165.

[2] Haji-Sheich A., Lakchimanarayan R., Lou David Y.S. Contined flow in a partially-filled rotating horizontal culinder // Trans. ASME: J. Fluid Eng. 1985. Vol. 106, N 5. P. 270-278.

[3] Шрагер Г.Р., Штоколова М.Н. , Якутенок В.А. и др. Моделирование вязкого течения со свободной поверхностью внутри вращающегося горизонтального цилиндра / / Теоретические основы химической технологии. 2005. Т. 39, № 3. С. 303-309.

[4] Johnson R.E. Steady-state coating flows inside a rotating horizontal cylinder //J. Fluid Mech. 1988. Vol. 190. P. 312-342.

[5] Fomin S., Watterson J., Raghunathan S., Harkin-Jones E. The run-off condition and hydraulic jump in rimming flow of a non-Newtonian fluid // Proc. of the ASME 2001 Fluids Eng. Division Summer Meeting. New Orleans, 2001. Paper FEDSM 2001-18186.

[6] Ross A.B., Wilson S.K., Duffy B.R. Thin-film flow of a visciplastic material round a large horizontal stationary or rotating cylinder //J. Fluid Mech. 1991. Vol. 39. P. 137-187.

[7] Fomin S., Hashida T., Watterson J., Fundamentals of steady-state non-Newtonian rimming flow // J. of Non-Newtonian Fluid Mech. 2003. Vol. 11. P. 19-40.

[8] Шрагер Г.Р., Козловродов А.Н., Якутенок В.А. Моделирование гидродинамических процессов при переработке полимерных материалов. Томск: Изд-во Том. ун-та, 1999. 219 с.

[9] Бреввиа К., Вроувел Л., Теллес Ж. Методы граничных элементов. М.: Мир, 1987. 524 с.

[10] Баттерфилд Р., Бенерджи П. Методы граничных элементов в прикладных науках: Пер. с англ. М.: Мир, 1984. 494 с.

[11] Шульман З.П., БАЙков В.П. Реодинамика и теплообмен в пленочных течениях. Минск: Наука и техника, 1979. 296 с.

[12] ЛАДЫЖЕНСКАЯ О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 288 с.

Поступила в редакцию 25 мая 2006 г.

Правила для Авторов

<http://www.ict.nsc.ru/mathpub/comp-tech/>

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

2. Все файлы предоставляются на дискете 3.5" формата 1440 Кбайт. Предпочтительнее пересылка файлов по электронной почте [email protected] в виде *.zip архива.

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

4. Статья должна сопровождаться разрешением на опубликование от учреждения, в котором выполнена данная работа. В сопроводительном письме необходимо указать почтовый адрес, телефоны, e-mail автора, с которым будет проводиться переписка.

5. Для каждого автора должна быть представлена (на русском и английском языках) в виде отдельного файла следующая информация:

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

о телефоны с кодом города (дом. и служебный), факс, e-mail, URL домашней страницы о область научных интересов (краткое резюме)

6. Материалы следует направлять по адресу: редакция журнала "Вычислительные технологии", Институт вычислительных технологий СО РАН, просп. Акад. Лаврентьева, 6, 630090, Новосибирск, 90, Россия, Игорю Алексеевичу Пестунову (отв. секретарь) — тел.: +7(383)3308785, e-mail: [email protected]; Галине Григорьевне Митиной (зав. РИО).

Рекомендации по оформлению статьи в LTEX

В редакцию следует направлять исходный файл, подготовленный в формате LTEX 2е в классе jctart (допускается использование стандартного класса article).

Файл класса jctart.cls можно скачать с сайта ЖВТ: http://www.ict.nsc.ru/win/mathpub/comp-tech/.

1. Структура файла в формате LTEX 2£:

\documentclass{jctart} \usepackage{amsmath}

\begin{document}

\setcounter{page}{1}\pagestyle{myheadings}

\markboth{<^0. Фамилия автора(ов)>}{<КРАТКОЕ НАЗВАНИЕ СТАТЬИ (ДО 40 СИМВОЛОВ)>} \title{<HA3BAHHE СТАТЬИ>\footnote{<Ссылка на поддержку (факультативно)>.}} ^uthor^sc^^^. Фамилия первого автора>}\\ \it{<MecTO работы первого автора>}\\ e-mail: \^{<Адрес первого автора>}\\[2тт] ^^^.О. Фамилия второго автора>}\\

\^{<Место работы второго автора (отличное от первого)>}\\ ...}

\date{}

\maketitle

\begin{abstract}

<Текст аннотации>

\end{abstract}

<Текст статьи>

\begin{thebibliography}{9}

<Библиография (\bibitem-список)>

{\small

\bibitem{} {\sc Иванов~И.И., Иванова~И.И.} К вопросу о вычислительных технологиях //

Вычисл. технологии. 1999. Т.~11, №-~11. С.~1123-1135. ...}

\end{thebibliography} \end{document}

(В конце файла даются:

<Перевод названия статьи на английский язык (или на русский, если статья на английском)> <аннотации на английский язык (или на русский, если статья на английском)>)

2. Список литературы составляется по ходу упоминания работы в тексте и оформляется по образцу:

Книга

Шокин Ю.И. Метод дифференциального приближения. Новосибирск: Наука, 1979. 222 с.

Бренстед А. Введение в теорию выпуклых многогранников: Пер. с англ. М.: Мир, 1988.

Рояк М.Э., СоловЕйчик Ю.Г., Шурина Э.П. Сеточные методы решения краевых задач математической физики: Учеб. пособие. Новосибирск: Изд-во НГТУ, 1998.

Finlaysön B.A. The Method of Weighted Resuduals and Variational Principles. N.Y.: Acad. Press, 1972.

Книга четырех авторов

Проблемы вычислительной математики / А.Ф. Воеводин, В.В. Остапенко, В.В. Пивоваров, С.М. Шур-гин. Новосибирск: Изд-во СО РАН, 1995.

Статья из продолжающегося тематического сборника

Федорова А.А., Черных Г.Г. О численном моделировании струйных течений вязкой несжимаемой жидкости // Моделирование в механике: Сб. науч. тр. / РАН. Сиб. отд-ние. Вычисл. центр. Ин-т теор. и прикл. механики. 1992. Т. 6 (23). С. 129-140.

Статья из журнала

Игнатьев Н.А. Выбор минимальной конфигурации нейронных сетей // Вычисл. технологии. 2001. Т. 6, № 1. С. 23-28.

Venkatakrishnan V. Newton solution of inviscid and viscous problems // AIAA J. 1989. Vol. 27, N 7. P. 285-291.

Труды конференции

Ivanov I.I. Problems in computational techologies // Intern. Conf. Comput. Techs. Novosibirsk, 1988. P. 225-229.

Препринт

Гуськов А.Е., Федотов А.М., Молородов Ю.И. Информационная система"Конференции". Новосибирск, 2003 (Препр. РАН. Сиб. отд-ние. ИВТ. № 1-03).

Диссертация

Деменков А.Г. Численное моделирование турбулентных следов в однородной жидкости: Дис. ... канд. физ.-мат. наук. Новосибирск, 1997. 123 с.

3. Иллюстрации вставляются в текст статьи с помощью команды \includegraphics, например:

\begin{figure}[htbp] \centering

\includegraphics{fig1.eps} \caption{<Подрисуночная подпись.>} \end{figure}

Наиболее предпочтительной формой представления иллюстраций являются файлы рисунков в векторном формате PostScript (.eps) или черно-белых растровых в форматах .bmp, .pcx, .tif с разрешением 300 dpi.

Все надписи на рисунках (обозначение осей и т.п.) должны быть выполнены в том же начертании (гарнитура "Roman"), что и в тексте статьи. Латинские символы — курсивом, из математической моды (x[fcj, z х 10-3, ф, P,...), цифровые обозначения на графиках — наклонно (№ кривой — 1, 2,...), цифры по осям — прямо (10, 15,...), единицы измерения — по-русски (кг, м,...).

Instructions für Authors

<http://www.ict.nsc.ru/math.pub/comp-tech/>

1. Papers may be submitted to the editorial board as two copies of the manuscript typed on one side of the standard A4 sheet (297x210 mm) and files of the manuscript in LTEX 2e format and files of the figures on a diskette.

2. All files should be submitted on a 3.5"floppy disc (1440 Kbytes) or sent to [email protected] as a *.zip - archive.

3. A separate page should contain a title, names of the authors, an abstract (not more than 300 characters) and keywords.

4. The paper should be accompanied by the publication permission from the organization, where the work was done. The enclosed letter should contain the postal address, phone numbers and e-mail of the corresponding author.

5. A separate file should contain the following information on each author:

o First name, second name, last name o Affiliation, position o Academic degree and title o Postal address

o Office and home phone numbers (including area code), fax number, e-mail address, homepage URL o Scientific interests (brief curriculum vitae)

6. All materials should be sent to the following address: Dr. Igor A. Pestunov (executive secretary), Journal of Computational Technologies, Institute of Computational Technologies SB RAS, Academician Lavrentyev Ave. 6, Novosibirsk, 630090, Russia. Phone +7(383)3308785, E-mail: [email protected]; Galina G. Mitina (publishing department manager).

Recommendations on submitting paper in LTEX

The source file should be submitted in LTEX 2e format using jctart class file (standard article class can also be used).

The files of appropriate jctart.cls class file can be downloaded from JCT web site: http://www.ict.nsc.ru/win/mathpub/comp-tech/.

1. The file structure in LTeX 2e format:

\documentclass[english]{jctart} \usepackage{amsmath}

\begin{document}

\setcounter{page}{1}\pagestyle{myheadings}

\markboth{<Name(s) of author(s)>}{<SHORT TITLE (LESS THAN 40 CHARACTERS)>} \title{<TITLE OF PAPER>\footnote{<Reference to supporting organization (optional)>.}} \author{\sc{<Name of the first author>}\\

\it{<Affiliation of the first author>}\\ e-mail: \tt{<Address of the first author>}\\[2mm] \sc{<Name of the second author>}\\ \it{<Affiliation of the second author>}\\ ...}

\date{} \maketitle \begin{abstract} <Abstract> \end{abstract} <Text of paper> \begin{thebibliography}{9} <References (\bibitem-list)> {\small

\bibitem{} {\sc Ivanov~I.I., Ivanova~I.I.} On computational technologies //

Computational technologies. 1999. Vol.~11, No.~11. P.~1123-1135. ...}

\end{thebibliography} \end{document}

2. A list of the references should be sorted according to the order of citations in the text and it should be written as in the following example:

Book

Finlayson B.A. The method of weighted residuals and variational principles. N.Y.: Acad. Press, 1972. Book by four authors

Problems of computational mathematics / A.F. Voevodin, V.V. Ostapenko, V.V. Pivovarov, S.M. Shurgin. Novosibirsk: SB RAS Publishing House, 1995.

Paper from continued subject transactions

Fedorova A.A., Chernykh G.G. On numerical modelling of viscous incompressible jet fluid flows // Modelling in mechamics: Scientific transactions / RAS. Siberian branch. Computing Center. Institute of Theoretical and Applied Mechanics. 1992. Vol. 6 (23). P. 129-140. Paper from journal

Venkatakrishnan V. Newton solution of inviscid and viscous problems // AIAA J. 1989. Vol. 27, N 7. P. 285-291.

Conference proceedings

Ivanov I.I. Problems in computational techologies // Intern. Conf. Comput. Techs. Novosibirsk, 1988. P. 225-229.

Dissertation

Demenkov A.G. Numerical modelling of turbulent wakes in homogeneous fluid: Dissertation for degree of candidate of physical and mathematical sciences. Novosibirsk, 1997. 123 p.

3. Figures should be included into the text using command \includegraphics{<figure file name>}, for example:

\begin{figure}[htbp] \centering

\includegraphics{fig1.eps} \caption{<Figure caption.>} \end{figure}

The preferred presentation form for illustrations is a figure file in vector format PostScript (.eps) or black and white bitmap formats .pcx, .bmp, .tif with 300 dpi resolution.

All figure inscriptions (axes definitions, etc.) should be done by the same font as in the text of paper ("Roman" type family). Latin characters should be done in italics in mathematical mode (x[k], z x 10-3, P,...), figures on axes — by straight font.

In papers, which are written in Russian, the units of measurement should be written in Russian.

В ближайших номерах/Forthcoming papers

Wang H., Cao D., Li S. Smoothing iterative method for a class of nonsmooth equations

Влнг Х., Дао Д. Ли С. Сглаживающий итерационный метод для одного класса негладких уравнений

Блохин А.М., Ибрагимова А.С. Об одном варианте метода прямых для симметрических t-гиперболических систем

Blokhin A.M., Ibragimova A.S. About one version a method of lines for symmetric t-hyperbolic systems

ЖумАгулов Б.Т., Темирбеков Н.М., Абдолдина Ф.Н. Определение степени загрязнения атмосферного воздуха города и моделирование процесса рассеяния вредных примесей

Zhumagulov B.T., Temirbekov N.M., Abdoldina F.N. Definition of a degree of pollution of atmospheric air of city and modelling of process of dispersion of harmful impurity

ЗоткЕвич А.А., ЛАЕвский Ю.М. Численное моделирование распространения ламинарного пламени на основе двухуровневых явных разностных схем Zotkevich A.A., Laevsky Y.M. Numerical modeling of laminar flame propagation based on two-level explicit difference schemes

Левин В.А., Луценко Н.А. Численное моделирование двумерных нестационарных течений газа через пористые тепловыделяющие элементы Levin V.A., Lutsenko N.A. Numerical modeling of two-dimensional time-dependent gas flow through porous heat-evolutional elements

Михеев С.Е. Метод точных релаксаций Miheev S.E. Method of exact relaxations

Паасонен В.И. Параллельный алгоритм построения гиперболических сплайнов

Paasonen V.I. Parallel algorithm for the constructing of hyperbolic splines

Пестунов А.И. Теоретические свойства теста "стопка книг" Pestunov A.I. Theoretical features of the "bookstack" test

Пузанов М.В., Шурина Э.П., Эпов М.И. Моделирование нестационарного электромагнитного поля методом векторных конечных элементов с использованием декомпозиции области

Puzanov M.V., Shurina E.P., Epov M.I. Modeling of time-dependent electromagnetic field by vector finite element method using domain decomposition

СтружАнов В.В. Об одном итерационном методе расчета напряжений в неод-носвязных телах

V.V. Struzhanov On one iteration method of stress calculation in non-simply connected solids