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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФИЗИКО-ХИМИЧЕСКИХ ПРОЦЕССОВ

УДК 519.6:532.62

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

АЛЬЕС М.Ю.

Институт механики УрО РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34

АННОТАЦИЯ. Исследуются эволюционные схемы численного решения краевых задач неньютоновских ползущих течений.

КЛЮЧЕВЫЕ СЛОВА: неньютоновские ползущие течения, краевые задачи, метод конечных элементов, эволюционные численные схемы.

ВВЕДЕНИЕ

Решение краевых задач неньютоновских ползущих течений сводится к решению систем нелинейных алгебраических уравнений большого порядка. Распространенными методами в таких случаях являются методы <асимптотической сходимости> [1], основанные на расширении исходной граничной задачи и сведении её к задаче Коши. Центральным моментом, при этом, является выбор подходящего линейного итерационного оператора перехода. Яркими примерами таких методов являются методы установления в газовой динамике [2, 3] и методы переменных параметров упругости и упругих решений в теории пластичности [4 - 6]. Воспользуемся идеями этого подхода для построения всегда сходящегося метода решения нелинейных краевых задач ползущих течений.

ПОСТАНОВКА ЗАДАЧИ

Решение задачи расчета медленного течения высоковязкой среды с движущейся свободной поверхностью, заполняющей в начальный момент времени t0 область V(/0),

сводится к определению в заданный момент времени t области V(V), скорости VI(ха,t),

давления р( ха, t), удовлетворяющих в V ^) уравнениям

-Ур + р& = 0, (1)

V и = 0, (2)

е- = +У и (3)

т = , (4)

при граничных условиях на поверхностях £ = и Бр П Бр =0

^ щ=и>;, (5)

Бр: щ (т} - р^ )пу = рь к, (Т - )щ = 0 (6)

и начальном распределении ц (t0). Здесь т}-, е}}. - тензоры вязких напряжений и скоростей деформаций; /лэ - эффективная вязкость; р - плотность; щ, ki - единичные векторы, нормальный и касательный к свободной поверхности; р0 - внешнее давление; gjj - метрический тензор; * - заданные распределения.

Свободная поверхность Бр движется в соответствии с кинематическим условием

д /

--= 0, где / - функция, описывающая свободную поверхность. Течения рассмат-

д1

риваются при следующих глубоких неравенствах для чисел Струхала Бк >> 1 и Рейнольдса Re << 1, капиллярного числа Са >> 1 - нестационарными и конвективными членами в уравнениях движения, силами поверхностного натяжения на свободных поверхностях, эффектом линий трехфазного контакта пренебрегается. На твердых стенках БЦ" и входных

границах 8ЦХ (З0 = ЗЦ" и 8ЦХ) задаются условия прилипания и профиль скорости. Зависимость эффективной вязкости описывается моделью Шульмана [7]

/ = (Т + (/У™) V, еи = (2е/е; ^, (7)

где т0 - предел текучести; / - динамический коэффициент структурной (пластической) вязкости; п, т - показатели нелинейности. Для сквозного расчёта при наличии квазитвердых областей течения (когда 0при т0 Ф 0), используется модификация [8]: е0 в (7)

заменяется на (еи + е) ( е ~ 10-4 .. ,10-3).

ПРОЕКЦИОННО-СЕТОЧНАЯ АППРОКСИМАЦИЯ

Для области V, ограниченной поверхностью З = Бр запишем слабую форму Галеркина для уравнений Стокса (1)

Ц/(еДУи + VV) -р^1 "\VyJV -\pgvJV - | /'уЖ = 0, (8)

V V Бр

и уравнения неразрывности (2)

I (Уи = 0 , (9)

V

где £ - компоненты вектора внешних поверхностных сил.

Произведя триангуляцию р и представив элементные распределения скорости и давления в виде (у - базисная функция) oi = уеои , р = ур1, 1еГп, из (8), (9) получим систему нелинейных алгебраических уравнений относительно неизвестных узловых значений скорости от и давления рт (т = 1, М )

N

Е^т I / (¿„хи/уу,+иуу) -у р, ^ ] у^п -

п=1 V

1 п

N N

-Ео\pgyrdV,,-^к от I /у/кя=0, (10)

п=1

N

Еот|(иУу= 0 , (т = 1,М; £,уё Гп, вёВп), (11)

п=1 V

^ п

где N, М - число конечных элементов и узлов в триангуляции р ; Гп - множество узлов элемента; Ь - множество элементов, у которых хотя бы одна грань используется при аппроксимации Бр; Вп - множество узлов элемента, находящихся на Бр; 0.ут - логическая процедура

1, если узел у е Гп конечного элемента п

О =

т

11, Уп е Ь

соответствует узлу т = 1,М триангуляции р ЛпЬ = <

I 0, в противоположном случае.

0 в противоположном случае;

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

Граничные условия на 8р (8) являются естественными для уравнений Стокса. Граничные условия на поверхностях удовлетворяются непосредственно путём

соответствующей коррекции матрицы и вектора правой части.

РЕШЕНИЕ СЕТОЧНОГО АНАЛОГА УРАВНЕНИЯ НЕРАЗРЫВНОСТИ И РАСЧЕТ ДАВЛЕНИЯ

Система уравнений (10), (11) является слабо связанной относительно давления. Воспользуемся методом искусственной сжимаемости [9], заменив (2) эволюционным соотношением

д^ + аУр' = 0, ^ = 0 р = 0 ц = 0, (12)

где д(-) - частная производная по (•); ^ - фиктивное безразмерное время; a > 0 - параметр

алгоритма, размерность которого совпадает с размерностью динамического коэффициента вязкости. Решение производится совместно с уравнениями движения и соответствующими граничными условиями. Физический смысл имеет только стационарное решение. Аппроксимируя (14) по неявной схеме с шагом А^ > 0 будем иметь

N Г ¿+1 я+1 1 N я

^1 ¥( р + аА^ и Уг¥( \¥( р1¥^п = 0 , (т = 1, М; £,ус Гп). (13)

п=1 V

Г

т

п=1 V

Параметры А<^, а подбираются численно, основываясь на оценке [9]: 4А^/(Де-1 + 0,5А^а)/Н2 < 1, где И - характерный шаг сетки. Уравнения (13) решаются итерационно совместно с уравнениями (10), в которых под узловыми значениями скорости и и давления р1 следует понимать их (я +1) приближения. Система за счёт ¡лъ (в0) является нелинейной.

ОСОБЕННОСТИ РЕШЕНИЯ СИСТЕМ НЕЛИНЕЙНЫХ КОНЕЧНО-ЭЛЕМЕНТНЫХ УРАВНЕНИЙ НЕНЬЮТОНОВСКОЙ ГИДРОДИНАМИКИ

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

Среди итерационных методов [10, 11] наибольшее распространение на практике получили одношаговые стационарные методы, которые, если представить уравнения (10), (11) в форме

И(б)б-G = 0, бТ = {и,р} , (14)

можно записать следующим образом [11, 12]

я+1 я я я

б=б-Ь-1(И (б) б- О). (15)

Основным требованием, предъявляемым к Ь, является локальная сходимость итерации (15) [12]. Если существует дбФ(б*) = дФ(б*)/дб (б* - точное решение системы (14)), то для локальной сходимости достаточно (и в существенном необходимо), чтобы спектральный радиус д матрицы удовлетворял неравенству

д(Е - Ь-1дбФ(б*)) < 1 . (16)

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

я я

Ь = д5(И (б) б) (17)

является квадратичная сходимость в е - окрестности точного решения [13, 14]. Трудность

0

состоит в выборе начального приближения 5 ее . Применение теорем о точке притяжения

0

[13, 14] для подбора 5 требует знания величин, вычисление которых для системы уравнений большого порядка является задачей, соизмеримой по трудности с исходной.

Исходя из ниже следующей оценки допустимого промаха е для скалярной функции

/(х) = 0 (|/(х)| >М1, |/(х}| <М2) х е (х*-е,х* +е) [15]

е < 2М1М 2-1

для элементных значений е0 для модели Шульмана имеем

е <

2е*

п 1 - п

1 + — + -

1/п

(18)

(19)

т

т Т" + (реи)

1/т

Из (19) видно, что е - окрестность, в зависимости от реологических параметров модели и реализуемого течения, может быть очень узкой, например, при т Ф 0 , т < 1 , п > 1 . В условиях пространственно-временной неоднородности решения, получение сходящегося процесса методом Ньютона связано со значительной вычислительной отработкой.

Более простым в практической реализации является метод простой итерации,

следующий из (15) при Ь = Н(5). В этом случае, помимо невырожденности матрицы Ь и требований к гладкости вектор-функции Ф(5), требуется выполнение дополнительных

условий, следующих из последовательности векторов ошибок щ = 5-5* [12]. С точностью до величин порядка 0(|Щ2) имеем

\Н-\5*)д3Н(5*)5*|| < 1 . (20)

При нарушении (20) сходимости нет, и даже при сколько угодно точном начальном приближении процесс простых итераций может разойтись. Для элементных значений е0 из условия (20) для реологической модели Шульмана имеем

пп

1--+ -

1/п

т

т т1п + /£)

1/т

< 1.

(21)

Из (21) следует, что метод простой итерации сходится только при определенных

*

сочетаниях реологических параметров модели и ео в конечных элементах:

если т = п- метод сходится при любых /, т0 (физически допустимых) если т Ф п,т0 = 0- метод сходится при п/т < 2

п

т/п

если т Ф п,т0 Ф 0 - метод сходится при — > 2, еи <

т

2

/л ^пт -2^

(22)

В окрестности точки п = 2т , при п / т ^ 2+ и фиксированном т0 Ф 0 , допустимое значение е*и возрастает и в пределе характер течения на область сходимости метода не влияет. При фиксированном п / т = 2+ и т0 ^ 0 допустимое значение е*и уменьшается, стремясь к нулю - метод сходится только для безградиентных течений. Однако, когда т0 = 0 влияние е*и на область сходимости метода исчезает, но решение сходится уже только при п / т < 2 . В наиболее практически важном случае т Ф п, т0 Ф 0 на параметры модели и реализуемые тензор-градиенты скорости накладываются очень жёсткие ограничения и метод не работает.

0

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

Расширим исходную задачу (1) - (6) заменив реологическую модель (4) следующим эволюционным уравнением

т, = 2аА/у + >0, ^ = 0, е, = 0, (23)

где ^ - фиктивное безразмерное время; ^ - параметр алгоритма, размерность которого совпадает с размерностью динамического коэффициента вязкости. Физический смысл имеет только стационарное решение (23). Аппроксимируя (23) по явной схеме с шагом А£/ > 0 ,

запишем следующий итерационный алгоритм (¡/ = ¡/^ А^-1):

я+1

У, тТ-Уг Р + р£ = 0,

я+1

я+1 1 я+1 я+1

= ¿(У и,+У, и), Уг и= 0,

я+1 я+1 я я я

т, = 2^// ( е , - е, ) + 2^3 (еи ) е,,

я+1 я+1

(Т - р g1J )n

я+1

= Г, и

= и.

(24)

Обобщённым решением задачи (1) - (6) назовём функцию VI , для которой справедливы уравнения (2) - (4), граничные условия (5) и которая удовлетворяет интегральному тождеству

Т (и)е, (р)ау = ау + \ ги ая (25)

V у Ър

для всякой достаточно гладкой функции ц, удовлетворяющей соотношениям (2), (3) и

однородным граничным условиям на поверхностях .

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

(и, и) = \ 2^е, (и)е, (б)ау, (26)

и нормой

является гильбертовым.

I I и I 1 = (и, и )1/2,

(27)

Теорема 1. Пусть существует единственное обобщенное решение задачи (1) - (6) параметр ¡Г алгоритма (24) таков, что выполняются неравенства

0 < т/^с2) <-!- д(^)е У) сис, < т^с2)

У

//

де

0 < т1 < т2 < 2,

(28) (29)

где с, - произвольный симметричный тензор второго ранга; 11 (с) = с - первый инвариант.

я

Тогда последовательность иг, определяемая системой (24), сходится к обобщенному решению задачи (1) - (6) со скоростью геометрической прогрессии

где

я+1 я

и - и

т = тах { 1 -

0

-

< т

{ I1 -т^,|1 -т2\}.

(30)

(31)

Доказательство, Эквивалентный итерационному процессу (24) алгоритм для обобщенного решения и. имеет вид

. £+1 Л . л . л

¡т"(и)е {и)йУ = \р£ и^У + | Г и. dS , (32)

где

s+1 s +1 s s s

T ( U ) = 2U f (e ( U ) - ev (u)) + 2U (e„ (U)K (u) •

- ] V - / -^ГГУ > ~гр

Вычитая из (32) аналогичное рекуррентное соотношение, записанное для приближения S, будем иметь

J 2U

s+1 s

e (и) - e (и)

ег] (u)dV = jj2U#

s s-1

e (и) - e (и)

s s

s-1 s-1

2^3 (e»^ (и) - 2^ (e»)ey (и)

^ (u)dV • (33)

£+1 X

Обозначим /лъ = /лъ(еи(и)). Принимая в качестве вектора и разность ц-ц , с учётом (30), (31), получим

s+1 s

и-и

E-ы -

s s s-1 s-1

u e (и) -U e (и)

s s-1

Uff ekl (и) - ek£ (и) 1

s s-1 s+1 s

2^ ^ (ekt (и) - ew (u ))(e, (и) - e„ (u))dV, (34)

где

Ekl = g'kg]t+gag]k )•

(35)

Воспользовавшись применительно к правой части (34) неравенством Коши-Буняковского будем иметь

2 Г У/2 / \1/2

[ ш £ X £-1 X X-1 1 [ ш £+1 £+1 £ £+1 £

<1 ¡л1]к1 2^еи(и- и)е"(и- и^У I ¡Л1]к1 2^еы(и-и)е"(о-о^У

s+1 s

и-и

s s

s+1 „ 1 д(Рз eii (и)) гДе Akt = Eki--

'ff

de

kl

Отсюда, с учётом Eijkle ev = I1(e ) и принятого условия (28), получим

, где m = max{|1 -mj,|1 -m2|}

s+1 s 2 s+1 s s s-1

и-и < m и-и и- и

(36)

При т1 + т2 < 2 имеем |1 ->|1 -т2| , а при т1 + т2 > 2: |1 -щ\<|1 -т2| . Отсюда, при т2 < 2 выполняется условие т < 1 и из (36) следует, что итерационный процесс (24) сходится к обобщенному решению ц задачи (1) - (6) как геометрическая прогрессия со знаменателем т

s+1 s

и-и

< m

s+1

0

и-и

(37)

Наименьшее значение величина m достигает при m1 + m2 = 2. Область допустимых значений параметра uff согласно (28), (29) имеет вид

Uff > 2 dev(u юо.

(38)

Теорема доказана.

Во введённом выше понятии обобщенного решения отсутствует гидростатическое давление. Это обусловлено тем, что функции и., ц. заранее (по определению) предполагаются удовлетворяющими уравнению неразрывности. Снимем последнее ограничение и введём обобщенное решение задачи (1) - (6) как вектор-функцию Ж = {и., р},

S

p

Л

Л

s

s

2

J

для компонент которой справедливы соотношения для тензора скоростей деформаций, граничные условия на поверхностях 5и = Б^1 и , определяющие уравнения и которая удовлетворяет интегральному тождеству

(Т (и) - р8г) е (и) - (Ур')р

(У = и гйУ +| Г игйл

(39)

для всякой достаточно гладкой вектор-функции Ж = { ц, р}, удовлетворяющей соотношениям (3) и однородным граничным условиям на 5 и . Из тождества (39) следуют уравнения движения, неразрывности и граничные условия на поверхности Бр . Определим скалярное произведение

Л

(Ж ) = |

(2/ (и) - р8м У (и) - (У и' ) р

йУ

(40)

и норму \\W\l = (Ж)1/2 , и покажем, что сходимость предложенного алгоритма (26) в норме | | и| | обеспечивает сходимость в норме ||Ж||.

Теорема 2. Пусть выполнены условия теоремы 1. Тогда последовательность ц, р определяемая системой (24), сходится в норме ||Ж|| к обобщенному решению задачи (1) - (6).

Доказательство. Эквивалентный итерационному процессу (24) алгоритм для обобщенного решения Ж имеет вид

s+1 s+1 Л ^ Л

П ( и , р К ( и) - (У и ') р

йУ = игйУ +| /' о^Б

(41)

s+1 s+1 s+1 s+1 s+1 s+1 s s s

где Пу ( и , р ) =Т ( и ) - р 8ц , Т ( и ) = 2 Л// (еу ( и ) - е ( и)) + 2 Лз (еи ( и))еу ( и) .

Вычитая из (41) аналогичное рекуррентное соотношение, записанное для приближения 5, получим

s+1

s s+1 s

2ц (еу ( и ) - еу (и) - ( р - р) 8<г

s+1

ег (и ) - ё ( и )

s s

еч (и ) - (Уг и' -Уг и' ) р(У =

Л

е^ ( и ) йУ.

s—1 s—1

2 ц е (и) - 2 ц е ( и )

Л s+1 s Л s+1 s

Принимая в качестве и = иг - ц., р = р - р, с учётом (40) и оценок для правой части тождества (33), будем иметь

2

s+1 s s s—1 s+1 s s s-1

(42)

Л+1 Л 2 Л+1 Л Л л-1 Л+1 Л 2 Л Л-1

Ж - Ж < т и - и и - и или Ж - Ж < т и - и

откуда и следует сходимость итерационного процесса (26) в норме Теорема доказана.

Полученные результаты справедливы для модели (4). Оценим область допустимых значений (38) параметра Л/ для модели Щульмана (7). Имеем

п

Иг > 0,5^з (еи)—

т

1 -

Тп + (ц (е„ )е )1/т

(43)

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

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

Л

Л

Л

Л

I

Б

Л

Л

Л' Л'

I

Б

р

Л

Л А

I

Л-1

Л

п

N Г

№ J

n=1 v

(T )ff — ^ Pi g

N Г

V1WJdVn — ^Qf J

n=1 V

T) ff — T

VPrdVn —

NN

—ZQm JPg'tyydVn — ZЛLQf J f!WßdSn = 0,

n=1

N

Y

m

n=1 V

v n

N

s+1

ZQf J (uX^ )wrdVn = 0,

n=1 V

(45)

где (Т)г = »(и} Уу +и; Уу), Т = »(еДи/ Уу + и'е Уу).

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

N

ZQm J

n=1 Vn

s+1

s+1

(tT ) ff +A £fauavavt g

f^^t aY it

vtvydVn — ZQf J

N

m

n=1 Vn

s s

) ff

(T) f — T+ y Pi g

N

r

m

n=1 Vn

—ZQf J Pg^rdVn — 2Л L Q f J f'WßdSn = 0.

^V iWydVn —

(46)

Давление на (£ +1) -итерации отыскивается из системы (13).

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

СПИСОК ЛИТЕРАТУРЫ

s

s

s

s

s

s

s

s

s

s

S

1. Марчук Г.И. Методы вычислительной математики. Новосибирск : Наука, СО АН СССР, 1973. 352 с.

2. Калиткин Н.Н. Численные методы. М. : Наука, 1978. 512 с.

3. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск : Наука, 1981. 304 с.

4. Ильюшин А.А. Пластичность. Основы общей математической теории. М. : Изд-во АН СССР, 1963. 272 с.

5. Биргер И.А. Некоторые общие методы решения задач теории пластичности // Прикладная математика и механика. 1951. Т. 15, вып. 6. С. 97-106.

6. Ворович И.И., Красовский Ю.П. О методе упругих решений // Доклады АН СССР. 1959. Т. 126, № 4. С. 740-743.

7. Шульман З.П. Конвективный тепломассоперенос реологически сложных жидкостей. М. : Энергия, 1975. 352 с.

8. Шрагер Г.Р., Якутенок В.А. Гидродинамика слива неньютоновской жидкости из осесимметричных емкостей // Инженерно-физический журнал. 1988. Т. 55, № 1. С. 66-73.

9. Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л. : Гидрометеоиздат, 1986. 352 с.

10. Самарский А.А., Гулин А.В. Численные методы. М. : Наука, 1989. 432 с.

11. Оден Дж. Конечные элементы в нелинейной механике сплошных сред. М. : Мир, 1976. 464 с.

12. Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М. : Мир, 1975. 560 с.

13. Канторович Л.В. О методе Ньютона // Труды Матем. ин-та им. Стеклова. 1949. Т. 28. С. 104-114.

14. Канторович Л.В. Некоторые дальнейшие приложения метода Ньютона // Вестник ЛГУ. Математика Механика. Астрономия. 1957. Т. 7, № 2. С. 68-103.

15. Самарский А.А. Введение в численные методы. М. : Наука, 1987. 288 с.

EVOLUTIONARY SCHEMES OF NUMERICAL SOLUTION IN NONLINEAR BOUNDARY PROBLEMS OF NON-NEWTONIAN CREEPING FLOWS

Alies M.Yu.

Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia

SUMMARY. The evolutionary schemes of computational solution in boundary value problems of non-Newtonian creeping flows are analyzed.

KEYWORDS: non-Newtonian creeping flows, boundary value problems, finite element method, evolutionary schemes of computational solution.

Альес Михаил Юрьевич, доктор физико-математических наук, профессор, главный научный сотрудник ИМ Уро РАН, тел. 8-912-856-38-24, е-таИ:аИв£ту@таИ.ги