Содержание к диссертации
Введение
Глава 1. Постановка задачи и предварительные сведения 9
1.1 Постановка задачи 9
1.2 Некоторые сведения из выпуклого анализа 10
1.3 Определение внешних аппроксимаций V, a, j и (Ро) . 13
Глава 2. Задача о течении среды Бингама в канале 16
2.1 Задача о течении вязкопластическои среды в канале . 18
2.2 Внешние аппроксимации вариационного неравенства для задачи о течении в канале 21
2.2.1 Первая внешняя аппроксимация HQ (П) 24
2.2.2 Вторая внешняя аппроксимация HQ(1) 25
2.2.3 Аппроксимация a,j,f 30
2.2.4 Сходимость внешних аппроксимаций 32
2.3 Сходимость ALG2 для задачи о течении в канале 34
2.4 Разностные схемы 43
2.4.1 Схема на разнесенных сетках 43
2.4.2 Схема на неразнесённых сетках 48
2.5 Численные результаты 49
2.6 Качественный анализ течения в трубах 55
Глава 3. Двумерная задача 59
3.1 Вариационная постановка для плоской задачи 59
3.2 Внешние аппроксимации для плоской задачи 60
3.2.1 Первая аппроксимация V 61
3.2.2 Вторая аппроксимация V 67
3.2.3 Аппроксимация aj,f 72
3.3 Алгоритм Узавы 76
3.3.1 Сходимость алгоритма 76
3.3.2 Первая разностная схема 81
3.3.3 Вторая разностная схема 85
3.4 Алгоритм ALG2 87
3.4.1 Сходимость алгоритма 88
3.4.2 Реализация алгоритма на разностных схемах . 93
3.5 Численные эксперименты 95
Глава 4. Трёхмерная задача 106
4.1 Внешние аппроксимации для пространственной задачи. 106
4.1.1 Первая аппроксимация V 107
4.1.2 Вторая аппроксимация V 109
4.2 Ядро дискретного оператора градиента 113
4.2.1 Структура дискретного оператора градиента 113
4.2.2 Тензорная структура оператора градиента и его ядра 124
4.2.3 Процедура ортогонализации 126
4.3 Спектр дополнения по Шуру 127
4.4 Стабилизация 131
4.5 Численные эксперименты 144
4.5.1 Аналитические решения для вязкой жидкости 144
4.5.2 Аналитические решения для вязкопластической среды 147
Глава 5. Моделирование нестационарных течений среды Бингама 149
5.1 Течение в канале 149
5.2 Течение в каверне 159
Заключение 160
Литература 164
- Некоторые сведения из выпуклого анализа
- Сходимость внешних аппроксимаций
- Первая разностная схема
- Ядро дискретного оператора градиента
Введение к работе
Актуальность работы
Существует широкий круг материалов, которые обладают поведением среды Бингама, а именно: ниже определённого предельного значения напряжений среда ведёт себя как жёсткое тело, выше этого предела - как несжимаемая вязкая жидкость. Примерами служат геоматериалы (высо-копарафинистая нефть, глины, грязи, сели, кристаллизующаяся лава), кровь в капиллярах, множество косметических и пищевых продуктов, строительные (свежий бетон, масляные краски) и химические материалы. При течении среды Бингама могут образовываться два типа зон: жёсткие зоны, в которых среда неподвижна или движется как твёрдое тело, и зоны течения. При этом естественным образом возникает «предельная» поверхность, разделяющая области с разным движением материала. Таким образом, характерной особенностью задач о течении вязкопластическои среды является необходимость строить решения в областях с неизвестной границей. Это обстоятельство создает значительные трудности при построении эффективных методов решения. Основная сложность при численном моделировании течения вязкопластическои среды связана с сингулярностью определяющих соотношений и невозможностью определения напряжений в жёстких зонах.
Существуют две основные группы методов для преодоления отмеченных математических трудностей: методы эффективной вязкости и методы, основанные на вариационной постановке. Методы эффективной вязкости (регуляризованные модели) состоят в аппроксимации недифференци-руемых определяющих соотношений гладкой функцией. Таким образом, среда рассматривается как нелинейная вязкая жидкость. Это упрощённый инженерный подход, который иногда приводит к нефизичным результатам, особенно для нестационарных задач.
Альтернативой регуляризованным моделям может служить вариационная формулировка, впервые предложенная для вязкопластическои среды А.А. Ильюшиным (1938). Одно из первых систематических исследований модели Бингама было предпринято П.П. Мосоловым и В.П. Мяснико-вым (1965-1967). Дюво и Лионсом (1972) с помощью выпуклого анализа задача минимизации функционала энергии сформулирована в виде вариационного неравенства, эквивалентного исходной системе уравнений.
Вариационная формулировка задачи автоматически учитывает наличие неизвестной границы, отделяющей области течения от жёстких областей. Более 30 лет назад в работах французских математиков 1 были предложены алгоритмы численного решения вариационных неравенств,
ХР. Гловински, Ж.Л. Лионе, Р. Тремольер (1976), М. Фортен (1984), П. Ле Таллек (1989)
однако они казались слишком сложными.В 2001 году алгоритм ALG2 был записан в виде удобной вычислительной процедуры, и стал активно применяться. Дискретизация проводилась методом конечных элементов, для которых он был предложен и обоснован.
Метод конечных разностей в рамках данного подхода не применялся, поэтому построение конечно-разностных схем, для которых обоснована аппроксимация и сходимость, является актуальной задачей, решению которой и посвящена диссертация. Необходимость разработки надёжных (обоснованных) и эффективных численных методов, соответствующих точной математической модели среды Бингама, обуславливает актуальность данной работы.
Цель работы
Основная цель работы состоит в разработке и исследовании новых численных методов решения задач вязкопластичности на основе вариационных неравенств с использованием разностных схем.
Методы исследования
Теория внешних аппроксимаций вариационных неравенств, теория разностных схем, методы оптимизации, методы вычислительной линейной алгебры.
Научная новизна работы
Представленные в диссертации результаты являются новыми и состоят в следующем:
Предложены две разностные схемы для смешанной постановки в переменных «скорость-тензор скоростей деформаций-тензор напряжений», являющиеся обобщением хорошо известных в вычислительной гидродинамике схем на разнесённых (МАС-схема) и полуразнесённых сетках, для задачи о течении в канале, плоской и трехмерной задач вязкопластичности. На основе теории внешних аппроксимаций обоснована сходимость решения конечномерной задачи к решению непрерывной.
Разработаны, обоснованы и численно реализованы новые методы решения задачи течения среды Бингама (на основе алгоритма Узавы и ALG2).
Для схемы на полуразнесённых сетках в трёхмерном случае проведено аналитическое исследование ядра дискретного оператора градиента. Предложены, обоснованы и численно реализованы два подхода
решения вырожденной задачи Стокса: поиск нормального решения на подпространстве, ортогональном ядру, и стабилизация разностной схемы путём введения дополнительного слагаемого в уравнение несжимаемости.
На основе разработанных разностных схем и известной схемы по времени проведено численное моделирование нестационарных течений вязкопластической среды для задач течения в каверне и каналах.
Практическая значимость работы
Предложенные разностные схемы могут быть использованы как для расчёта течений среды Бингама, так и других неньютоновских жидкостей, в том числе моделей с нелинейной вязкостью.
Апробация работы
Результаты диссертации обсуждались на научных семинарах Института вычислительной математики РАН, («Вычислительные и информационные технологии в математике» рук. проф. В.И. Лебедев, проф. Ю.М. Нечепуренко, чл.-корр. Е.Е. Тыртышников, 2008, 2009), Института математического моделирования РАН (2008), Института автоматизации проектирования РАН (рук. акад. О.М. Белоцерковский, 2008), Институт прикладной математики РАН им. М.В. Келдыша (семинар им. К.И. Бабенко, 2009), Института проблем механики РАН им. А.Ю. Ишлинского (рук. акад Д.М. Климов и акад. В.Ф. Журавлёв, 2009), на научно-исследовательских семинарах кафедр механико-математического факультета МГУ им. М.В. Ломоносова: теории упругости (рук. проф. И.А. Кийко, 2008), механики композитов (рук. проф. Б.Е. Победря, 2008), вычислительной математики (рук. проф. Г.М. Кобельков, 2008), оптимального управления (рук. проф. В.М. Тихомиров, 2009), теории пластичности (рук. чл.-корр. Е.И. Ломакин и акад. И.Г. Горячева, 2009); на семинаре кафедры высшей математики физического факультета МГУ им. М.В. Ломоносова (рук. проф. В.Ф. Бутузов, 2009), семинарах университетов Humboldt University, IWTM Kaiserslautern (2007), University of Cape Town, TU Cape Town, TU Darmstadt, Hausdorff Institute of Mathematics (HIM, Bonn) (2008).
Результаты диссертации докладывались на Международной конференции «Ломоносов» (Москва, 2005, 2009), «Ломоносовские чтения» (Москва, 2006, 2009), III International conference «Computational methods in applied mathematics CMAM-3» (Минск, 2007), II International conference «Matrix methods and operator equations» (Москва, 2007), Седьмой Всероссийский семинар «Сеточные методы для краевых задач и приложения» (Казань, 2007), Workshop «Viscoplasticity: from Theory to Application» (Тичино, Швейцария, 2007), International Conference on Numerical Analysis and Applied
Mathematics ICNAAM (2008), ENUMATH (Уппсала, Швеция, 2009), Workshop on Advanced Computer Simulations for Junior Scientists (Санкт-Петербург, 2009), Workshop «Viscoplasticity: from Theory to Application» (Лимассол, Кипр, 2009).
Публикации
По теме диссертации опубликовано 8 статей в рецензируемых журналах, из них [2-8] - в журналах из списка, рекомендованного ВАК.
Структура и объем работы
Некоторые сведения из выпуклого анализа
Пусть X — линейное пространство. Множество Л с X называется выпуклым, если для любых двух точек хт и Х2 из А и любого числа 0 из отрезка [0, 1] элемент 0хі + (1 — 0)х2 принадлежит А. Функцией будем называть отображение f : X — М, где R = Ж U {—со, +оо} - расширенная числовая прямая. С каждой функцией можно связать два множества: }, называемые соответственно эффективным мнооюеством и надгра-фиком функции f. Функция f называется выпуклой, если ее над-график - выпуклое множество. Функция f называется собственной, если f(x) —оо Vx, f(x) ф- +оо тождественно (domf 7 0)- Собственная функция выпукла тогда и только тогда, когда она удовлетворяет неравенству Иенсена: Строгая выпукостъ функции означает, что f ((1 — 9)xi + 9х2) (1 — 9)f(x!)+9f(x2) Vxbx2 ЄХ, 9 Є [0,1]. Функция f : X —» R называется полунепрерывной снизу на X, если она удовлетворяет одному из двух равносильных условий: Vex є R {х є X I f (x) x} замкнуто, Vx0 Є X liminfx_ Xof(x) f(x0). Пусть V - рефлексивное банахово пространство (с нормой ), К - непустое замкнутое выпуклое подмножество V и функция J : К —- Ж, J- выпуклая, собственная, полунепрерывная снизу. (1.2.2) Рассмотрим задачу минимизации infj(u). (1.2.3) Любой элемент v Є А такой, что J(v) = inff(u). (1.2.4) называется решением задачи (1.2.3). Теорема 1 Пусть функция ] - выпуклая, собственная, полунепрерывная снизу. Если J коэрцитивна на К, а именно limj(v) = +00 при и — оо, и Є А, (1.2.5) или множество К ограничено, (1.2.6) то задача, (1.2.4) имеет хотя бы одно решение. Если J строго выпукла на
К, то решение единственно. Рассмотрим следующую задачу минимизации: infj(u). (1.2.7) Предположим, что J(u) можно представить как верхнюю грань по р некоторой функции [и, р), то есть (Z - рефлексивное банахово пространство) J(u) = sup {и, р) Vu є V. (1.2.8) peZ Поскольку всякая выпуклая собственная полунепрерывная снизу функция совпадает с верхней гранью семейства всех не превосходящих ее аффинных непрерывных функций, можно, вообще говоря, записывать J (и) в виде (1.2.8), причем часто несколькими различными способами. Тогда задача (1.2.7) примет вид: inf sup(u, р). (1.2.9) Двойственной к (1.2.7) естественно назвать задачу sup inf (u, р). (1.2.10) рє2 ueV Пусть функция L определена на множестве А х В и принимает конечные значения. Пока можно считать, что А и В - два произвольных множества. Предложение 1 Если функция L определена на множестве А х В и принимает конечные значения, то sup inf L(u, р) inf sup L(u, p) (1.2.11) Пара (v, q) называется седловой точкой функции на А х В, если (v, р) (v, q) {и} q) Vu Є Л Vp є В. (1.2.12) Предложение 2 Функция, определенная на А х В и принимающая конечные значения, имеет седловую точку (v, q) на Ах В тогда и только тогда, когда 2 (v, Ч) = max inf (u, р) = minsup(u, р). (1.2.13) рєВ иєА иєА VEB 2Отметим, что замена inf (sup) на min(max) означает, что нижняя (верхняя) грань достигается. На пространстве V рассмотрим функцию J(v) следующего вида: T(v) - la(v, v) + j(v) - (f, v). (1.2.14) Здесь a(u, v) - некоторая билинейная непрерывная форма u, v — a(u, v); из непрерывности следует, что a(u, v) cuv, с = const, u, v є V. Форма a(u, v) предполагается симметричной, следовательно, a(u, v) = a(v, u) Vu, v є V. (1.2.15) Функционал v — j(v)- выпуклый, собственный, полунепрерывный снизу из V — R, (1.2.16) (f, v)- скалярное произведение элементов f Є V и v є V, где V -пространство, сопряженное к V. Для того чтобы J(v) удовлетворяла (1.2.5), достаточно выполнения следующего условия (v-эллиптичность): существует a 0, такое, что a(u, u) au2 Vu Є V. (1.2.17) Задача infj(v), veV. (1.2.18) эквивалентна вариационному неравенству найти v є V, такое, что a(v, u-v)+j(u)-j(v) (f, u-v) VueV. (1.2.19) 1.3. Определение внешних аппроксимаций V, a, j и (P0) Практическая реализация численных алгоритмов предполагает сведение задач к конечномерным. Для этой цели существуют два достаточно общих класса методов: методы внутренней аппроксимации и методы внешней аппроксимации.
Теория внутренних аппроксимаций применяется, как правило, при исследовании схем МКЭ или каких-либо других методов типа метода Галеркина. Основным признаком внутренней аппроксимации является принадлежность аппроксимирующего пространства V исходному V, что позволяет исследовать близость приближенного и точного решений в V. При использовании конечно-разностных схем аппроксимации исходных задач (уравнений или неравенств) приближенным решением является сеточная функция. Удобным и достаточно простым способом восполнения является построение кусочно-постоянной функции VH и кусочно-постоянных восполнений разностных производных. Однако в этом случае VH не принадлежит V и теория внутренних аппроксимаций не может быть использована. Здесь используется теория внешних аппроксимаций. Суть ее состоит в том, что рассматривается более широкое, чем V, пространство F и сравнивается некий образ crv элемента v в пространстве F с неким образом рьУк элемента VH В F. Систематическая теория внешних аппроксимаций впервые была развита Сеа [12], частный случай аппроксимации «крест» был введен Ж.-Л. Лионсом [36], исследование для случая уравнений проведено в работе Тремо-льера [53].
Исследование внешних аппроксимаций для задачи Стокса приведено в [110]. Аппроксимация V Зададим гильбертово пространство F и оператор сг, обладающий свойствами о" Є (V, F), о - инъективный оператор из V в F. (1.3.1) Рассмотрим далее семейство конечномерных пространств Ун, снабженных (гильбертовыми) нормами инук. Предположим, что заданы операторы продолжения -ръ. Є (VH, F). Операторы продолжения называются устойчивыми, если их нормы ограничены раномерно по Тг: IIPHIU(VH.F) С. (1.3.2) Семейство VH реализует внешнюю аппроксимацию пространства V (внешняя аппроксимация называется сходящейся), если ин Є VH, PhiiH - E, слабо в F =Ф В, є aV (1.3.3) и если Vu Є V, Зин Є VH, такое, что PHU-H —» сої сильно в F и UHVH С. 1 (1.3.4) В качестве ин часто берут гін = гни, где r : V — VH - оператор сужения. Аппроксимация а Зададим теперь билинейную форму CLH(UH, VH) на VK и предположим, что Ыин, vh) cuhYhKvh, (1.3.5)
Сходимость внешних аппроксимаций
Теорема 5 Пусть v - решение задачи (Рон)к к = 1,2; v -решение задачи (Ро); если К — 0, то V Эх" Эх" f силъпо 6 (L2(n))3» (2.2.51) о -) ) ( 3v 9v 1 ч {qHv, Tiqhv, нЧ) - 4 v, —, — силъко в {ІГ[Ж))Л. (2.2.52) Доказательство. В этом разделе будем опускать индекс k = 1,2 для двух различных аппроксимаций. В формуле (2.2.12) положим uh = 0: fh (VK) dx; cih (VH, VK) + jh (VH) отсюда, а также из неравенства jh. (VH) 0 И (2.2.14) следует оценка IKIIvk = СІК (Vh, Vh) fL2(Q)qHVhL2(Q) f L2[Q) vhvh V Tl, следовательно, IKIk C VH. (2.2.53) Учитывая неравенства (2.2.14),(2.2.53) llq Vhlk n) С, Vh, 5iqhVhL2(Q) C, Vh(i = 1,2). (2.2.54) В силу ]Ph.U(vH,F) С (2.2.17),(2.2.23) и предыдущих неравенств имеем IIPHVKIIF С. (2.2.55) Из последовательности Vh можно извлечь подпоследовательность, по-прежнему обозначаемую VH , такую, что PHVH — В, слабо в F, (2.2.56) что в совокупности с (2.2.18),(2.2.24) показывает, что L, = crw, следовательно, PhVh. — o v слабо в F. (2.2.57) Но в силу (2.2.45),(2.2.46),(2.2.47) из (2.2.57) вытекает, что liminf ah(vh, vh) a(w, w), (2.2.58) liminf jh(vh) j(w). (2.2.59) Выберем произвольное u Є HQ(Q) , тогда, согласно (2.2.19),(2.2.25) Б uh (. = Thu) є Vh, такое, что PhUh — сої сильно в F. Кроме того, из условия PhUh. — cat сильно в F и в силу (2.2.57) из условий (2.2.43), (2.2.44) следует, что ah(vh, Uh) - CL(W, U), (2.2.60) из (2.2.48) jh(uh)- j(u), (2.2.61) из (2.2.50) - что (fh, UH-vh) - (f, u-w). (2.2.62) Перепишем (2.2.12) ah (VH, vH) + jh (Vh) CLH (vh, uh) + jh (uh) - (fh, uh - vH). (2.2.63) Используя в неравенстве (2.2.63) полученные результаты (2.2.58) -(2.2.62), находим a(w, w) + j (w) a(w, u) +jn(u) - (f, uh —w) Vu є Hj(O), следовательно, w = v, что позволяет утверждать, что т ъУъ. — crv слабо в F. Для доказательства (2.2.51), (2.2.52) рассмотрим элемент VH Є VH, для которого Ри%і — o"vv сильно в F, (2.2.64) и положим Хн = ah (VH - vhj vh - vh). Имеем Хн = ah (vh, vH)-aH (VK) Vh)-ah (vH vh) + aH (vh, vh). Из (2.2.63) следует, что Хн ан (vh, йк) + jh (uh) - jh (vh) - (fh, uH - vH) -- aH (vh, VH) - ah (vh, vH) + aH (vh, vh) Используем теперь (2.2.43), (2.2.44) и получим lim sup Хн a (v, u) + j(u) — J (v) — f (u — v) — a (v, v). Полагая теперь u = v, получим, что Хн —» 0, откуда, ввиду (2.2.42), имеем vh-VHvh - 0, однако в силу phU(vh)F) С (2.2.17),(2.2.23) PhVh — PhVh — 0 сильно в F, откуда и следует (2.2.51), (2.2.52). Напомним, что Vh = { vh I vK = (vij)MeaH, v Є R}. Рассмотрим первую аппроксимацию. Определим пространство Qh = (чн qh = (Чн) Чн) clh = { ЧЧ+ІА МЄЩ» % = ( чу+ідімєп } (2.3.1) Дискретный оператор градиента VIH : VH —» Qh определяется как VIHUH = (6IUH, 6214,.) Таким образом, пространству QH принадлежат Qh, Uh, (VihUh).
Рассмотрим вторую аппроксимацию. Определим пространство Qh = {qn qn = (qk ) qk = ( qW.j+i/iWns, 1 = 1,. (2.3.2) Дискретный оператор градиента Vjh Мг — QH определяется как V2HUH = (f 1 uh, TT2UH) Пространству Q принадлежат qh, цК) (V2uh) Введем дискретные аналогиx функционала J (u, q) и расширенного лагранжиана т: Jh(uh, qh) = Ы2 + qh2 + (Tsih(qh) - (fh, uh)(2.3.3) rh(uh) qH M-H) = JH(UH, qn) + (Uh, VKuh - qh) + - Vhu,v -- qh{2.3.4) с г 0. Параметр т не входил в исходный функционал J(UH) , а в расширенный лагранжиан входит. Необходимо проверить, что тройка (v, у, Л), являющаяся седловой точкой лагранжиана rh., не зависит от выбора г. Другими словами, для произвольного неотрицательного г эта тройка также будет седловой точкой г/н. Теорема 6 Пусть (vh, yh, Ан) - седловая точка А-иСи-н, qh М-к) Тогда она также является седловой точкой г/н для любого г 0 и VH - -решение задачи минимизации функционала J (u). Доказательство. Пусть (VH, УН, AH) - седловая точка rh. на VH х QH х Qh. Тогда Vuh Є VH, qh Є Qh, М- Є Qh справедливо Th(vh, Th, Hh) AhK, У h лн) rh(uh, qh, Ah). (2.3.5) Докажем, что для седловой точки (VH, УН, AH) выполняется ун = VHVH- Рассмотрим первое неравенство (2.3.5): rh(vh, Ун, М-н) AH(VH, УН, АК). Так как в левую и правую части неравенства входят одинаковые величины VH и у н, то оно равносильно ( Hhi vH h - Тн) (Ah, VHVh - Ун) V(iH Є QH =Ф VhVh = уh. Следовательно, (Ah, VhVh - у J + т/2/д VhVh - ун2 dx = 0. Рассмотрим второе неравенство (2.3.5). Так как оно выполняется Vqh,Uh, возьмем произвольный элемент ин Є VH И ДЛЯ него выберем qh так, Дальнейшее изложение справедливо для обеих аппроксимаций, поэтому соответствующие индексы опущены. что qh = VHUH. Поэтому (Лн, VhUh-qh) +r/2jn VhUH-qh2 dx = 0. Тогда a/2vH2 + /2yH2 + JH(yH)-(f,VH)=/:rH(vH) yh, AH) rh(uH, qn, Ah) = ос/2ин2 + n/2qh2 + jh(qh) - (f, uh), следовательно, Jh(vh, yH) Jh(uh,qh) VuH Є Vh, qh = VhUh, поэтому Jh(vh) ЫШг) Vuh Є Vh- Таким образом, первая компонента седловой точки (vh, yh, Ah) лагранжиана Th.("u-h, qh М-н) является также точкой минимума для Jh(Uh). Рассмотрим теперь расширенный функционал r h для произвольного фиксированного г 0. Убедимся, что (vh, yh Ah) удовлетворяет неравенствам r h(vh, yh» Vh) r h(vh, yh, лк) A- h(-LLh, qh, Ан). (2.3.6) Так как yh = VhVh, то r h(vh, yh, Цн) = Л- н іг, yh Ah) V h є Qh. Первое неравенство (2.3.6) доказано, перейдём к доказательству второго неравенства (2.3.6). Возьмем произвольный элемент Wh Є Vh и рассмотрим Uh = tWh+ (1 )vh, где 0 t 1 (очевидно, что Uh Є Vh), выберем qh = yh = VhVh- Тогда правое неравенство (2.3.5) принимает вид x/2twh +(1- t)vH2 - a/2vh2 + t(f, vh - wH)+ +Tt2/2Vh(wh-Vh)2 + t(Ah, Vh(wh-Vh)) 0 Vwh Є Vh. Учитывая выпуклость функции x/2,2, Vwh Є Vh ( tWh +(1- t)vh2 tWh2 + (1 -1) IKII2, f t(K24vh2)+t(f, Vh-Wh)+ t2Vh(wh-Vh)2+t(Ah) Vh( h-Vh)) 0. Разделим на t и устремим t — 0: j (IKII2-IKII2)+ (f, vh-Wh) + (Ah, VhK-Vh)) 0 VwhG VK. (2.3.7) Аналогично, выберем u = VH и возьмем произвольный элемент TJ н є Qh- Положим qh = УнП- "t)+trih, 0 t 1, и из второго неравенства (2.3.5) получим ІІУн(1- +іт,н2 + )н(ун(1- +іЛк)- ун2-Зн(ун)+
Первая разностная схема
Через UH, VH будем обозначать пространства сеточных векторных функций, заданных только во внутренних точках Пі и П2, Vh = Uh х V . Введём дискретное пространство для тензорных функций: QIK = (qn I qn = Ык , Чн » Чн . 4 (q V )tj = ЧК(ХІІ), (чн2)у = qi?( ij) ч e a3, (qi2)y = qJi2Uii). (q?)y = й1 ), - є a4} Определим разностный аналог оператора дивергенции VIH- : Vh — Рн и тензора скоростей деформаций Diu : VK —) Q1K: (VIH vK)ij = (uy - иі-і,,0Лч + (vy - vi b1 )/h2, (Dih(Vh))i,j = (Uy-Ui-ijVHb (D K(vH))y = (vg-vy_1)/h.2, (DjH(vH))i,j = (D2h(vH))y = (uij+i kj)/2H2 + (vi+1J - vy)/2h!. Для скалярных, векторных и тензорных сеточных функций зададим скалярное произведение формулами [Vh, ФтОїк = Pi,7 Pi,jHiH2 (wh, vh)ih = (X.wljvlj + Hw2.ivy)hlh-2 a, a2 (qn, AH)IH = (Ц(ЧЦАЦ + qgAg) + JjqgAg + q A2])) . Далее в тексте будут использоваться следующие нормы, все они вычислялись на основе введенных скалярных произведений: Рнт = (Рн,Рн)!н vhih = (vh,vH)!h2, AKih = (Лн, Ан)]н2-Рассмотрим дискретный аналог лагранжиана (и; ц ) (3.3.2): ГУ IH(UH) ЛН) = 2IWIm + HPih(uh)2h+o"s(лн. Dlh(uH))ih- (fh, uh)ih (3.3.26) Определим замкнутое выпуклое множество Лік пространства QIK A1h = {АніАк Є Qm, Ан = Ан,ЛнІ 1} норма тензорной сеточной функции: ІАнІч = \УЮ2 + (А22)2 + ( й/2,ж/2)2 + (ла,/2,)+./2)2+ + 1)2 + (А?;2)2 + (ла1АН1/2)2 + (HlV2i +V2Y-+ + уТл;,1 +(л?2) +№У1]_Ч2у-+(л2;1/2)+" + IXhlui Д)+і/2 = 7 Ід/ОД )2 + (А22)2 + (АЦ1/2,,+1/2)2 + (A2;,/2ij+1/2)2+ i+1/2,j+1/2 )2 + (AU,/2 ,)+1/2 + V(A«+.)2 + (A2j+,)2 + (A,/2,W/2)2 + (A2I,/2 (3.3.27) Норма, введенная таким образом, соответствует внешней аппроксимации функционала. Чтобы найти седловую точку дискретного лагранжиана к(и.н, Ah) применяется дискретная версия алгоритма (3.3.16)-(3.3.18) в следующей форме: Пусть задано произвольное Ан є Ацг; Для n = 0,1,2,... выполняется: ШАГ 1: Найдем у как решение системы (3.3.28) осчі - \±Amv + V1HPK = o-s VIH AS + fh, Vlh-vS = 0; ШАГ 2: Вычисляем (3.3.29) AH+1=PAH H + r DlhK)). Если A +1 — AJJ e, повторить с шага 1.
Таким образом, первый шаг алгоритма заключается в решении обобщенной задачи Стокса, второй представляет собой метод проекции градиента на выпуклое множество (сферу). Для введеной дискретизации с учетом VIH VH = 0 выполняется 2(D1h(vk),Dm(vh)) = (V1hvh) VIHVH) = (-AIHVH.VH), (3.3.30) при этом получим следующие сеточные аналоги оператора градиента VIH : РН - UH Х VH И дивергенции Vm- : Qm - UH х VK: (VihPh)i,j = ((РІ+1.) -Ру)Аь (py+i -Py)/H2), (3.3.31) ТП __П т12__12 т21__21 _22 _22 _ (V,h тн)у = ( + П ZbLT?=Li + Г ) - (3-3.32) оператору Лапласа Лін : VH — UH Х Ун соответствует стандартный пятиточечный шаблон «крест» на сетках JQI и 1 Первый шаг заключается в решении обобщенной задачи Стокса - с постоянной матрицей системы и переменной правой частью. Первыми итерационными методами решения стационарной задачи Стокса в переменных «скорость-давление» были методы Эрроу-Гурвица и Узавы. Несмотря на то, что с момента создания этих методов прошло более пятидесяти лет, они остаются основой для создания новых итерационных методов [119, 2]. Рассмотрим систему сеточных уравнений для стационарной задачи Стокса Бели система невырождена, то при решении методом Узавы сначала решается уравнение для давления SHPH = ф SnVh = VH Ай1 VHPH = Vh A FK = p, а затем вычисляется скорость VH = (Ан)_1 (FH — VHPH)- Условие LBB-устойчивости для разностной схемы задачи Стокса можно записать следующим образом. Предположим, что 0 с h-ih "1 С с некоторыми константами с и С. Тогда существует константа Со 0, не зависящая от h = max{Hi, h.2}, такая что sup VH-VH) Сорн11 VPHGPH. (3.3.34) VhGUhxVh (-AHVH,VHJ2 Оператор S)v = VH-A 1 VH является самосопряженным на Рн, неравенство. (3.3.34) эквивалентно Соїн $н, ГАе Ъг - единичный оператор на Рн-
Для полученной разностной схемы условие (3.3.34) выполняется [31]. = N 1 для заданных натуральных Ni,N2. Определим следующие сеточные области: П3 = {xtj = (irn, jh2) і = 0,..., Nb j = 0,..., N2}. Определим пространства компонент сеточных функций скорости и давления V2h =-К) = (Ui,j,Vy) = (U(xlj),v(xij)) I Xij Є П3, V0,j = Vi,N2 = VNl,j = vii0 = 0}, P2h = {ptj := p(xij) I xtj Є П4, Pij = 0} у Через V"2h будем обозначать пространства сеточных векторных функций, заданных только во внутренних точках С1-\ . Пространство сеточных тензорных функций обозначим через Q2H Qih = (qn I qn = {qi1, qi2, qj?, q? 1; (qhkj = qh(x4), хц є n2). Определим разностный аналог оператора дивергенции V2K- V —» Р2н и тензора скоростей деформации D2H : V2h — Q2H: (V2H vh)i,j = (ui+ij+i - uy+i- + iti+i,j - uy )/2Hi + + (Vi+ij+i -vi+1J+vy+1 -vy)/2H2, (3.3.35) (D2h(vh))y = (гн+ij+i -Uy+i +Ui+ij -Uy)/2Hb (D2J(vH))i,j = (Ui+i j+i - -щ+i.j + uy+i - -Щ,0/41г2 + + (Vi+ij+i -vy+1 +vi+1 j -vi(j)/41 . (3.3.36) Определим скалярное произведение и нормы: (рн, Phbh = "руФірНіНг, n3 (wh)vK)2H = Ц(І,ІЧІ + ЦіЧі Н2 Wh = (WH,WH)T, vH = (vh,v2)T, S, (qn, ЛнЬн = L№4] + ОЙ + ЧЙЛЯ + 03) 2, llPhlbh = (Ph,Ph)2(2, vh2h = (VK,VH)2(2, ЛнІ2к = (Ah, Лн)2н2 Дискретный лагранжиан имеет вид /"V 2h(Uh, Лн) = у IMllh +M-D2h(Uh)2h+"s( Л Н D2H(uh))2h - №г, Uh)2h. (3.3.37) с использованием введённых в этом разделе норм. Определим замкнутое выпуклое множество Лгк пространства Q2H А2Н = { НІ Ан Є Q2hi Ан= Ah, Ah. 1}, Ahi,j = ((AjJ )у+2(Л{ )?5Н-(А )? )1/2 - норма тензорной сеточной функции.
Для рассматриваемой дискретизации применение алгоритма (3.3.28)-(3.3.29) приводит к использованию следующих сеточных аналогов оператора градиента Vih Ргн — V2H И дивергенции V2H- : Q2H — V2H (V2hPh)i,j = ((РУ - Pi-u + Pij-i - Pi-i.j-i )/(2Hi), (pg -py-1 +Pi-i,j -Pi-i,j-i)/(2H2)), (3.3.38) (V2h-Ah)ij = ( 2h + 2b 2hi 2h-2 и оператора Лапласа Д2н У ь. 2н (2-2.29) случае hi = h2 = К для Аги получим так называемй «косой крест» (2.2.30). Хорошо известно, что для задачи Стокса построенная разностная схема с операторами (2.2.29),(4.4.2),(3.3.38) является LBB-неустойчивой. Чтобы убедиться в этом, достаточно заметить, что сеточный оператор градиента V2H (4.4.2), в отличие от непрерывного, имеет на пространстве Р2н нетривиальное, ядро вида Ker(V2h) = span(p\p2), pjj = (-1) -1, р, = (-1)і+і+1-1, (3.3.39) поэтому оператор SIK имеет нетривиальное ядро на Р2н и неравенство Coin $2К вьпіолнено только с Со = 0. Следовательно, конечно-разностное LBB условие устойчивости (3.3.34) не выполнено, однако этот недостаток вполне компенсируется вторым порядком аппроксимации и удобством использования разностной схемы. В случае, если Szh вырожденна, потребуем выполнения условия (p_!_kerV2h Тогда задача разрешима, одним из стандартных методов решения системы (3.3.33) является поиск нормального решения уравнения S2HPK = Ф на подпространстве пространства Р2н, ортогональном Ker(V2n) (3.3.39).
Ядро дискретного оператора градиента
Как было отмечено ранее (раздел 3.3.3 главы 3), схема на полуразнесённых сетках является LBB -неустойчивой уже в двумерном случае. Это не является препятствием для её применения, но, несомненно, создаёт определённые сложности при использовании данной схемы. В предыдущей главе был описан метод поиска нормального решения рассматриваемой системы, который может быть применён и в трёхмерном случае. Поскольку для реализации данного подхода должны быть известны базисные векторы kerBh, данный раздел посвящен изучению размерности и структуры ядра дискретного оператора градиента 4.2.1. Структура дискретного оператора градиента. Для нахождения размерности воспользуемся методом сингулярного разложения. SVD является наиболее надежным методом, позволяющим вычислить размерности и определить базисы образа и ядра линейного оператора [68]. Сингулярным разложением действительной та х п -матрицы А называется всякая ее факторизация вида где V - ортогональная та х та-матрица, W - ортогональная п х п-матрица,а Z -диагональная та х и-матрица с неотрицательными диагональными элементами. Величины О"І ЯВЛЯЮТСЯ сингулярными числами матрицы Л. Ранг произвольной матрицы Л равен рангу диагональной матрицы І в ее сингулярном разложении. С помощью SVD было проведено исследование ядра оператора В [78]. Для плоского случая размерность ядра была равна 2. Для трехмерного случая имеет место следующая зависимость: Тс = 3N — 1 , где N - число узлов по каждому из направлений по p(N — 2,... ,6), к. -размерность ядра В. Исследуем структуру ядра оператора В (разностный аналог оператора градиента) в плоском случае, то есть найдем ненулевые векторы Pij Є RNp, удовлетворяющие уравнению Вр = 0: Введем следующие обозначения: оператор Вх соответствует , Ву - -. Первое уравнение (4.2.1) соответствует Вхр = 0, второе - Вур = 0. Мы будем использовать обозначения В, Вх, Ву как для операторов, так и для соответствующих им матриц.
В - прямоугольная матрица размера Nu х Np. Каждая строка оператора В представляет собой запись разностного уравнения в точке с координатами (i,)). При этом в уравнение будут входить 4 соседние точки р с соответствующими координатами. Для удобства в дальнейшем будем рассматривать дополнительную сетку с узлами в точках, в которых задано значение р . Мы можем рассмотреть квадрат с вершинами в этих точках, причем точка с координатами (i,j) будет являться центром квадрата. Назовем его клеткой. Рассмотрим одну из клеток. Этой клетке соответствуют два уравнения (4.2.1) . Пусть в узлах клетки стоят числа Для того чтобы набор чисел, стоящих в узлах сетки, принадлежал ядру оператора В, должны удовлетворяться уравнения Эти равенства накладывают ограничения на числа a,b,c,d: a = d, b = с. To есть клетка будет иметь вид Вид должна иметь каждая клетка. Т.к. каждое ребро клетки (i, j) входит в соседние клетки, то, присоединяя соседние клетки, получим Таким образом, ядру оператора В принадлежат векторы вида (4.2.1), где a, b - произвольные числа. То есть размерность ядра равна 2 и базисом могут быть следующие векторы: Применение метода SVD подтвердило данное утверждение. Для матриц различной размерности и четности числа узлов размерность ядра была одинаковой и равна 2. Перейдем к трехмерному случаю. В пространственном случае узлу (i, j, к) соответствует три уравнения (для u, v, w). В каждое уравнение входит значение р в вершинах «кубика», центром которого является точка (г, j,k). Определить общий вид ядра для трехмерного случая намного сложнее. Можно попытаться построить несколько «частных», решений, воспользовавшись полученным решением для плоского случая. В качестве верхней грани кубика возьмем клетку Чередуя эти сетки в разной последовательности, будем получать элементы ядра. Элементы каждой из этих сеток удовлетворяют «плоским», уравнениям Вхр = 0 и Вур = 0 (4.2.2), поэтому первые два уравнения (4.2.7) выполняются, т.к. в каждое уравнение входят элементы двух соседних слоев, а для каждого слоя эти уравнения выполняются автоматически.
Третье уравнение будет выполняться, если сетки повторяются или чередуются(под единицей в верхнем слое стоит нуль в нижнем слое, и наоборот). Аналогично можно чередовать плоские сетки в направлениях yz и xz. Построим базис для числа узлов N = 23 (2 узла по каждому направлению).Мы должны найти значения в узлах. Для 8 неизвестных должны выполняться 3 уравнения (4.2.7). Ранг этой системы равен 3, следовательно, размерность ядра равна 5. Можно указать следующие "кубики": п В первых шести случаях плоские сетки повторяются по одному из направлений, в следующих двух - чередуются. Сумма векторов в каждой паре равна единичному вектору (е ?), поэтому в базис может входить только один вектор из пары. Это 4 вектора. В качестве пятого вектора достаточно взять второй вектор из любой пары или единичный вектор е ?. Например, можно взять векторы (еі, ез, є5, Є7, е8). Таким образом, мы нашли базис. SVD дает тот же результат (dim kerB = 5)
