CC BY
25
3/2006
ПОЛУИТЕРАЦИОННЫИ МНОГОСЕТОЧНЫЙ МЕТОД ДЛЯ РЕШЕНИЯ МНОГОМЕРНЫХ КРАЕВЫХ ЗАДАЧ РАСЧЕТА КОНСТРУКЦИЙ
А.Б. Золотов, М.В. Белый, В.Е. Булгаков, В.Г. Бельский
Введение. Численное решение трехмерных краевых задач теории упругости и теплопроводности, встречающихся в строительстве и технике, вызывает порой существенные трудности даже при использовании современных высокопроизводительных ЭВМ. Это объясняется тем обстоятельством, что доступные пользователю сеточные разбиения не всегда позволяют обеспечить требуемую точность аппроксимации при решении реальных задач со сложными контурами, краевыми условиями и т.д.
При использовании итерационных методов, например, основанных на экономичных разностных схемах, как правило, наблюдается плохая сходимость по низкочастотным составляющим решения, в особенности при решении задач в областях, сочетающих массивные и вытянутые части. В задачах теории упругости это приводит к невыполнению интегральных условий равновесия. Между тем известно, что для большинства уравнений эллиптического типа низкочастотные составляющие решения носят плавный характер и для их сеточной аппроксимации не требуется много узлов. Более того, каждому диапазону частот соответствуют составляющие решения определенной плавности и своя сетка. Поэтому кажется естественным построить итерационный процесс на последовательности сеток, самая мелкая из которых является исходной, а остальные - вспомогательными. Решение задачи на самой крупной сетке целесообразно осуществить прямым методом. При сравнительно небольшом числе узлов эффективен и двухсеточный подход.
Некоторый опыт применения вспомогательной крупной сетки был накоплен еще до широкого применения вычислительной техники, в частности при решении задач механики [13]. Интенсивное развитие многосеточные методы получили начиная с работ Р.П. Федоренко и А.Б. Золотова, как в нашей стране, так и за рубежом (см., например, [3,7-9,12,14]).
В настоящей работе строится эквивалентный оператор для итерационного процесса, рассмотренный ранее в [1,2,5]. Оператор сочетает последовательность сеточных разбиений и удобен для программных реализаций в виде универсальных программных комплексов.
Рассмотрим задачу
Аи = /, (0.1)
где и, / - сеточные функции, принадлежащие пространству Н со скалярным произведением ( • , • ); А - вариационно-разностный или конечноэлементный самосопряженный положительно определенный оператор (матрица).
Для решения задачи (0.1) будем использовать двухслойную итерационную схему
Вик+1 = Вик - тк (Лик - /), (0.2)
где В - эквивалентный оператор, который должен быть легко обратимым и обеспечивать хорошую сходимость (0.2); Тк - итерационный параметр на к-й итерации.
Если и - точное решение задачи (0.1), то для вектора погрешности гк, очевидно, имеет место соотношение
гк+1 = Ткгк, где Тк = Е-ткВ~1Л ; гк = ик - и. (0.3)
Идея многосеточного полуитерационного метода состоит в том, что для построения эквивалентного оператора В используются вспомогательные вложенные укрупняющиеся сеточные разбиения. Рассмотрим сначала двухсеточный вариант полуитерационного метода, а затем и общий случай.
1. Двухсеточный вариант полуитерационного метода.
1.1. Построение эквивалентного оператора. Наряду с исходной мелкой сеткой £ введем вспомогательную крупную сетку £0, являющуюся подсеткой сетки £.
Пространство сеточных функций на £0 со скалярным произведением ( • , • ) обозначим Н0. Элементам Н0 будем приписывать индекс 0, например, Ц)^Н0. Пусть Q - оператор интерполяции (продолжения) сеточных функций с крупной сетки на мелкую, а Q* - сопряженный с ним оператор сужения с мелкой сетки на крупную, удовлетворяющей условию
(и Оуй) = ((2*и Vо)о, и е Н, V е Н0 .
Определим сужение оператора Л на крупную сетку:
Л0 = Q*ЛQ . (1.1)
Эквивалентный оператор, определяющий процесс (0.2), строится следующим образом. Представим оператор Л в виде Л = Б + & , где Б - легко обратимый оператор, & = Л - Б . В качестве Б можно выбрать, например, верхнюю или нижнюю треугольную часть матрицы Л, включая главную диагональ, или Б = diag(а11, ..., апп), где ац - элементы главной диагонали матрицы Л.
Поставим задачу, «близкую» к исходной задаче (0.1) и определенную на исходной мелкой сетке:
БУ + GQvо = /, (1.2)
где У0 - решение следующей задачи на крупной сетке:
^ = и /0 = . (1.3)
Запишем близкую задачу в операторной форме:
Ву = / или у = В- /.
Выражая V из (1.2), (1.3), получаем
V = Б->(Е - GQЛ(-1Q*)f ,
откуда
В = П1(Е - GQЛ(-1Q*). (1.4)
Легко видеть, что обращение оператора В состоит из решения задачи (1.3) малой размерности (на крупной сетки) и выполнения ряда других нетрудоемких операций. Следовательно, этот оператор является легко обратимым.
Умножим (1.2) слева на Q* и, учитывая (1.3), разрешим полученное уравнение относительно Уо:
Уо = (А, - Q*GQУQ*Dv . (1.5)
Подставим (1.5) в (1.2)
^ + GQ(Aо - Q*GQ)D]v = /.
Отсюда получаем эквивалентный оператор в виде
В = [Е + GQ(Aо - Q*GQ) . (1.6)
Непосредственной проверкой устанавливается, что В1 AQ = Q, т.е. оператор В ~1А действует на векторы вида V = Qv0, У0 е Н0 как тождественный. Отсюда следует, что оператор В близок к А на множестве медленно изменяющихся, плавных сеточных функций, которые могут быть хорошо приближены интерполянтами с крупной сетки, что позволяет надеяться на быструю сходимость итерационного процесса (0.2) по низкочастотным составляющим решения.
1.2. Действия на одной итерации для случая двух сеток. Для наглядности подытожим алгоритм одной итерации:
„к+1
и = и - ТкУ , где ук = D-V - GQvk0); У0к = ; гк = Аик - /.
Вопрос поиска Тк будет рассмотрен ниже.
1.3. Исследование сходимости. Введем пространство сеточных функций на мелкой сетке НА со скалярным произведением [и, v] = (Au,v) и нормой
|1и11 =и], . Рассмотрим вопрос о сходимости итерационного процесса (0.2) по норме пространства НА.
Пусть итерационный параметр Тк = 1. Тогда оператор перехода (0.3) с учетом (1.4) имеет вид:
Тк = Е - В- А = (Е - D- А]Р ,
где Р = Е - QAlQ*А . Для любых и, v е НЛ справедливы равенства
(17)
(1.8)
[Р 2и, v] = [Ри, v]; [Ри, v] = [и, Pv].
Следовательно, Р - оператор проектирования. Отсюда ||Р||=1. Используя (1.1) и (1.8), можно показать, что
PQ = 0, &АР = 0.
Значит Р: НА ^ Жи ядро оператора Р имеет вид КегР=¥, где V- подпространство интерполянтов с крупной сетки, т.е. векторов типа v=Qvо, Vо е Но, а Ж - ортогональное дополнение к V в НА. Таким образом, второй сомножитель в (1.7) гасит плавные, низкочастотные составляющие вектора погрешности, тогда как первый сомножитель осуществляет погашение высокочастотных компонент.
Оценим норму оператора перехода. Используя свойства оператора Р, получаем
|(Е - D - А)и\
и ^ и || тки || |
|| тк ||= = виР-
иеНА || и || иеЖ
■<|| Е-D- А |
и |
(19)
А.Б. Золотое, М.В. Белый В.Е. Булгаков, В.Г. Вельский 3/200^^^^^^ГСуТНИК
В правой части неравенства (1.9) стоит норма оператора перехода двухслойного итерационного процесса с эквивалентным оператором, равным И. Отсюда следует, что при = 1 полуитерационный метод сходится не хуже, чем другие методы, определяемые видом оператора £), причем можно надеяться на лучшую сходимость полу итерационного метода, так как подпространство Ж уже, чем НА. Сходимость метода может быть улучшена за счет подбора величин итерационных параметров т/г Ниже в пункте 3 для одномерной модельной задачи будут установлены величина скорости сходимости и оптимальное значение параметра.
2. Полуитерационный метод в общем случае.
2.1. Построение эквивалентного оператора. Введем последовательность сеточных разбиений 8р,р =0, 1..., да, расположенных в порядке возрастания количества узлов, т.е. - самая крупная сетка, - самая мелкая сетка, совпадающая с исходной, причем £ ■[ для любогор =1, 2..., да является подсеткой сетки Л'/г Пространство сеточных функций на со скалярным произведением ( • , • )р обозначим Нр. Элементам Нр будем приписывать нижний индекс р, например, ир е Нр. Введем операторы интерполяции Ор. Нр_| —р =0, 1..., да, и сопряженные с ними операторы сужения Ор, удовлетворяющие условию
("р, орур-а )р = (ор"р, V1 > 11 р е нр > V е Нр ' •
Обозначим
ip =
E, если р = m,
Qm ■■■Qр+l, еСЛИ Р = 0, 1 m-1,
т.е. 1р - оператор интерполяции с сетки Sp на сетку Sm (самую мелкую).
Определим сужение оператора исходной задачи (0.1) на вспомогательные сетки:
Ар = IpAIp, р = 0, 1, ..., m . Здесь Am = А. Легко видеть, что для любого р =0, 1..., m-1
Ар = Qp+1 Ар+Qp+1 . (2.1)
Каждый из операторов Ар, р =0, 1..., m-1, представим в виде Ар = Dp+Gp, где Dp - легко обратимый оператор.
По аналогии с двухсеточным случаем поставим задачу, близкую к исходной и определенную на исходной сетке, с помощью цепочки равенств
А0У0 = /0 , DpVp + GpQpVp-1 = fp, р = 1, 2, ..., m , (2.2)
где /р = I* f. Вектор v = Vm получается из (2.2) последовательным определением Vp, Р =0, 1..., m:
v0 = А-/0 , vp = D~p fp - GpQpVp-1X Р = 1 2, m .
Отсюда получается рекуррентная формула для оператора В 1 = В)я1, обратного оператору близкой задачи:
(2.3)
в-1 = А-1
К = В7(Е? -), Р = 1, 2, ..., т.
Здесь Ер - тождественный оператор в пространстве Нр.
Умножим второе равенство в (2.3) на Qp и, учтя, что $/ = / -1 = В
выразим Ур_1'.
Подставим (2.4) в (2.3): Отсюда
V1 = (Вр-1 Вур.
[Вр + Ор$р (Вр-1 - )-%Вр у = /р .
'р-1 р-1'
(2.4)
В- = А
Вр = [Ер + Ор$р(Вр-1 -арОДр]Вр , р = 1, 2, ..., т.
(2.5)
Построенный оператор близкой задачи В = Вт будем использовать в качестве эквивалентного оператора итерационного процесса (0.2).
Покажем, что оператор В ЛА действует на интерполянты с самой крупной сетки ¿0 как тождественный, т.е.
В-1 А1- = /-. (2.6)
Действительно, легко проверить, что В11 А1$1 = .
Пусть В-Ур-$р-,..а = $р-1...$1 .
Тогда с учетом (2.1), (2.5) имеем:
В-1 АДр...$1 = В-АДр..& -Б-р1(Ар -Бр№рВр\А^^ = . (2.7)
Последовательно применяя (2.7), получаем (2.6), что и требовалось. Из (2.6) следует, что, как и в случае двух сеток, оператор В близок к А на множестве плавных сеточных функций, которые могут быть хорошо приближены интерполянта-ми с самой крупной сетки.
Действие на вектор оператором В-1 складывается из решения задачи на самой крупной сетке и поэтапного перехода на мелкую сетку с помощью нетрудоемких операций интерполяции, умножения на Вр и Ор. Задача на крупной сетке, как правило, может быть эффективно решена в классе точных методов (методов исключения). Следовательно, оператор В является легко обратимым.
2.2. Действия на одной итерации для случая т сеток. Вычисляется
где у* = В-1« -GpQpvkp-l), р =0, т;
„.к+1 „.к _ к и = и - Т,
ук = А-и *гк у0 " А0 /0Г
* к
К =1 /
гк = Аик - / .
Вопрос поиска т. будет рассмотрен ниже.
2.3. Исследование сходимости. Введем пространства НА сеточных функций на сетках ¿р, р =0, 1..., т со скалярными произведениями [ир,ур]р = (Аирур)р и нормами ||ир || = [ир, ир ]°р'5. Рассмотрим вопрос о сходимости итерационного процесса (0.2) по норме пространства Нт= Н^ (индекс для нормы в пространстве Н^ будем опускать).
Пусть итерационный параметр хк = 1. Тогда оператор перехода (0.3) с учетом (2.3) имеет вид:
Тк = Е - В-1 А = (Е - Б? А)(Е - ОВЖ А). Отсюда, используя тождество
Е - 1рВ-;1*рАр = (Е - 1рБ-;1;А)(Е - /р_,ВГ^Гр- А) ,
получаем
Тк = Хт ... Го, где Г0 = Е - ¡о А-1/о* А ; tp = Е - 1рП-;1РА, р =1, 2..., т. (2.8) В пространстве НА рассмотрим подпространства Ур, р =0, 1..., т, состоящие из векторов вида V = /рур, ур е Н^, т.е. Ур - подпространство интерполянтов с сетки 8р на исходную сетку. Очевидно, У0 е... е Ут = НА. Обозначим через Жр ор-
тогональное дополнение к Ур в НА.
В пункте 1 было показано, что Г0 - ортогональный проектор на Жд. Значит, оператор Г0 в (2.8) гасит интерполянты с самой крупной сетки.
Для всех р =1, 2..., т
Гр/Р = ¡р (Ер - Ар) . (2.9)
Следовательно, Ур является инвариантным подпространством для оператора Гр.
Для любого w е Жр имеет место равенство /*рА~№ = 0. Поэтому если w е Жр, то = w. Следовательно, Жр также является инвариантным подпространством для оператора Гр, причем на векторы из Жр этот оператор действует как тождественный.
Таким образом, операторы Гр,р =1, 2..., т отвечают за погашение тех компонент вектора погрешности, которые являются интерполянтами с соответствующих сеток В частности, оператор Гт отвечает за погашение самых быстроизме-няющихся составляющих погрешности.
Оценим норму оператора перехода:
т
НТкРПНГтН.
р=0
Так как Г0 - проектор, то || Г0|| = 1. Дляр =1, 2..., т-1
Г
|| гр ||= тах 1, Бир-
vеУn
| Л
|| V | )
Используя (2.9), получаем
||^|| || tp¡pVp || ||/р (Ер - о; Ар )Vp|| А
Бир и и = Бир и Т и = Бир —-—^—^^—— =|| Ер -Ор Ар ||р
II VII г^Н-А \\1рУр\\ грЕН$ ||/^р|| Поэтому
т-1
II Тк ||<|| Е-И-'А II Пт^х(1, II Ер-В-;Ар ||р ). (2.10)
р=1
Оператор /)р, р =1, 2..., т можно выбрать так, чтобы выполнялись условия || Ер ~0~рАр ||р<1, р = 1, 2, ..., т . (Это возможно, так как Ер ~Б~рАр - оператор перехода интерполяционного процесса с эквивалентным оператором равным 1)р для решения задачи на сетке Тогда || Тк\\ < 1, откуда следует, что интерполяционный процесс (0.2) сходится. 113
Заметим, что полученная оценка (2.10) является очень грубой и не может быть использована для вычисления величины скорости сходимости, которая на самом деле, как правило, оказывается существенно лучше, чем получающаяся из (2.10).
3. Анализ сходимости метода.
3.1. Введем на отрезке [0, 1] последовательность сеток ¿р, р =0, 1..., т с постоянными шагами Ир = 2т-рМ_1, где Мкратно 2т. В качестве исходной задачи на мелкой сетке рассмотрим систему вариационно-разностных уравнений | А_1[ - и(1 -1) + 2и(/) - и(I +1)] = /(I), I = 1, 2, ..., М-1 [ и (0) = и (М) = 0.
Операторы и 0р* определим следующим образом:
для любого I = 1, 2, ..., Мр -1
Г ир-1 (0.51), если I четное, ©А^ХО = | 0.5{ир-1[0.5(1 -1)] + ир-1[0.5(1 +1)]}, если I нечетное
для любого ] = 1, 2, ..., Мр-1 -1
($>р)(/) = 0.5[ир(2]-1) + 2ир(2/) + ир(2] +1)] .
Сужения исходной задачи на вспомогательные сетки ¿р, р =0, 1..., т, имеют
К[ - ир(I -1) + 2ир(I) - ир(I +1)] = /р(I), I = 1, 2, ..., Мр-1 ир (0) = ир (М) = 0.
Собственные значения и соответствующие им собственные векторы для оператора нар-й сетке имеют вид:
вид
Xkp = 4h-1 sin2 (0.5knM;), Фр(i) = 2-05 sin(kn/M-1), k, i = 1, 2, Mp -1. Используя тождества
QpФр-1 = cos2 (0.5knMp Vp+1 - sin2 (0.5knMp^MP-k ; (QP+I ... QÍФ"т)(i) = 4(hA)-1 sin2 (0.5knM-1)205 sin(kп/M-1) ,
можно показать, что в случае двух сеток, т.е. при т=1,
-»-»-i > k i * kп . 2 kп i k kп . 2 kп м-k
B 1 ^ =1 1 - cos—sin2-I mf + cos — sin2-фМ k
1 I M 2M J 1 M 2M 1
B-1 ^-k =-cos—cos2 kп ~k ■ 1 1 ■ 1 M
—Ф1+11+cos—cos2 — |m1-k 2M 1 I M 2M J 1
Отсюда следует, что линейная оболочка векторов ф 1 и ф 1/_к является для В _1А инвариантным подпространством, что позволяет легко найти собственные значения оператора В -1А, а следовательно, и собственные значения оператора перехода (0.3), которые в рассматриваемом случае имеет вид:
= 1 — т[1 + cos2 {ккК4 1)], k = \, 2, ..., М0 ,
причем кратность собственного значения |д=1-т равноМ0. Известно [11], что норма вектора погрешности асимптотически убывает быстрее геометрической прогрессии со знаменателем р = шах |. Наилучшая сходимость имеет место при
к
X = xopt = 2[2 + cos2 (ккМ'1 )]-1.
При этом
р = Ртш = cos2 (tzM~1 )[2 + cos2 (тiM-1 )]_1 <1/3 . В случае трех сеток (т=2) можно показать, что собственные значения оператора перехода = 1 - хЪ,к , где Ъ,к - корни уравнения
+16 sin2 (2а) cos2 (2а)[1 + cos2 (4а)](^ -1) - sin2 (8а)} =0, а = 0.5кШ~1.
Таблица 1
m M Pmin Topt
1 20 0.3279001 0.6720998
2 20 0.3280000 0.6719999
1, 2 40 0.3200000 0.6679999
1, 2 80 0.3330001 0.6669998
Вычисление спектра оператора перехода и определение ртщ в случае т=2 проводилось численно.
Как видно из таблицы 1, в случаях одной и двух вспомогательных сеток сходимость практически одинаковая, что говорит о предпочтительности многосеточного варианта в целях снижения размерности задачи на крупной сетке.
3.2. Для автоматизированного выбора величины итерационного параметра при решении задач по схеме (0.2) удобно использовать метод минимальных невязок [10], который приводит к соотношению
т, = (АВ У, гк)(АВ V, АВ У Г1 , где гк = Аик -/ .
Наряду с двухслойной итерационной схемой (0.2) применялась трехслойная схема
ик+х = ик -акВ 1 (Аик -/)-р,(ик -Ик '), (3.1)
где параметры и [3/. также выбираются из условия минимизации нормы невязки на шаге.
В ряде случаев сходимость итерационного процесса (3.1) может оказаться лучше, чем у (0.2). Проиллюстрируем сказанное на примере решения модельной плоской задачи теории упругости (плоская деформация £=1000, \'=0.3) с помощью авторской программной реализации метода. Расчетная область с приложенной внешней нагрузкой, условиями закрепления и нанесенной сеткой конечных элементов показана на рис. 1. Для построения эквивалентного оператора использовалась последовательность сеток (6x6 узлов), ^ (11x11 узлов), (21x21 узлов). Сетка выделена на рис. 1 жирными линиями. В качестве р = 1,2 ис- _ _ _
___115
пользовались диагональные матрицы, совпадающие с главными диагоналями матриц Ар, p = 1,2. В табл. 2 приведены данные о сходимости двухслойной и трехслойной схем для рассматриваемой задачи. В этой таблице и ниже используются следующие обозначения:
|| u || 1 = max| ui |; || u || 2 = max| ui |; || u || 3 = (u,u)0 5;
(\\Auk - f\\ t /\\Au0 - f\\ t) -100%, i = 1, 2, 3
4. Опыт авторских программных реализаций и решения реальных задач.
В связи с большим объемом вычислений при решении трехмерных задач важно построить эффективные алгоритмы, реализующие основные операции, встречающиеся как в итерационном процессе, так и при подготовке к нему. К таким операциям относятся умножение вектора на матрицу, сужение сеточных функций и матриц на более крупное сеточное разбиение, интерполяция сеточных функций на более мелкую сетку, формирование матрицы исходного дискретного оператора. При программной реализации указанных операций для современных ЭВМ оказалось эффективным организовывать выполнение основных арифметических операций и операций обменов с внешними запоминающими устройствами не над отдельными переменными и без организации циклов средствами языка высокого уровня, а над большими массивами чисел с помощью специальных программ, написанных на языке Assembler.
В разработанном авторами программном комплексе дискретизация краевых задач осуществляется по методу стандартной области [3,4,6]. Задача всегда решается на сетке, топологически эквивалентной прямоугольной с размерами М^М^хМз по трем индексным направлениям. Коэффициенты исходных дифференциальных уравнений имеют свои истинные значения на исходной области и продолжаются нулем вне ее. Такой подход позволяет также упорядочить поток исходной информации и сделать ее задание удобным для пользователя. Использование ряда приемов позволяет задавать пространственные области сложной формы. При нумерации неизвестных сначала изменяют индексы узла, а затем - номер неизвестного в точке, если у краевой задачи имеется более одного неизвестного в точке (задача теории упругости). Тогда матрица A и ее сужения будут иметь блочную структуру, причем все блоки A¡j одинакового размера и каждый из них содержит максимум 27 ненулевых диагоналей, размер центральной N = М^М^хМз,, а остальные близки к этому. С помощью таких диагоналей целесообразно хранить матрицу во внешней памяти.
Рисунок 1
Таблица 2
№ итерации Двухслойная схема Трехслойная схема
^ 2 83 Тк 82 83 а к вк
1 100.00 100.00 0.644 100.00 100.00 0.644 0.00
2 303.71 49.73 0.854 303.71 49.73 0.872 -0.0811
3 149.94 15.44 0.632 155.34 13.784 0.711 -0.00724
4 93.67 7.79 0.895 83.26 6.67 0.991 -0.150
5 54.55 4.37 0.796 41.61 3.49 0.950 -0.218
6 34.47 2.72 0.825 19.138 1.522 1.143 -0.102
7 22.61 1.780 0.842 9.334 0.778 0.962 -0.155
8 14.91 1.18 0.811 5.000 0.430 0.979 -0.169
9 10.21 0.825 0.897 3.126 0.281 1.264 -0.341
10 6.95 0.573 0.786 1.818 0.153 1.098 -0.177
11 4.97 0.416 0.967 1.016 0.087 1.105 -0.211
12 3.46 0.300 0.762 0.580 0.049 0.968 -0.156
13 2.54 0.224 1.046 0.356 0.0316 1.198 -0.283
14 1.81 0.166 0.737 0.216 0.0184 1.058 -0.209
15 1.35 0.126 1.123 0.127 0.0109 1.132 -0.201
16 0.999 0.0947 0.716 0.108 0.0079 0.822 -0.121
17 0.751 0.0723 1.187 0.0576 0.0045 1.118 -0.195
18 0.588 0.0544 0.700 0.0366 0.0028 1.266 -0.312
19 0.425 0.0415 1.244 0.0219 0.0017 1.139 -0.247
20 0.327 0.0313 0.687 0.012 0.0009 1.033 -0.208
21 0.256 0.0241 1.182 приближение к предельной точности, определяемой длиной четырехбайтового слова
22 0.199 0.0182 0.721
23 0.145 0.0140 1.180
24 0.114 0.0107 0.715
25 0.085 0.0082 1.197
Умножение вектора на матрицу осуществляется следующим образом. Пустьу=Аи, тогда
V- = ¿1, 1 = 1 2, .... I ,
1=1
где I - количество неизвестных в точке.
Через ак обозначим к-ю ненулевую диагональ блока Агу, представленную в виде вектора, дополненного нулями до размера N в начале, если к > 0, и в конце, если к < 0; при этом -13 < к < 13 и а° - центральная диагональ. Тогда
- 13
V- =111 з = 1 2 1 ,
1 =1 к=-13
где компоненты векторов Щ определяются следующим образом: при к > 0
К (я) =
при к < 0
К (я) = ■
аз (я)и1 (я + ак X 1 < я < N -ак 0, N > я > N -а к
0, 1 < Я < а к
а-к (Я -ак >] (Я), N > Я >ак.
Здесь щ - величина сдвига, зависящая от номера ненулевой диагонали, которую легко вычислить. Как правило, < N.
Алгоритмы умножения реализованы с помощью специальны« процедур на языке Assembler. Матрица в виде ненулевык диагоналей, исходный вектор и результат их умножения хранятся во внешней памяти ЭВМ. При этом в каждый момент в оперативной памяти должны располагаться минимум три массива размером N.
Для построения эффективного алгоритма интерполяции последняя производится не по пространственным переменным Xj, X2, X3, а по целочисленным ij, ¿2, /3, являющимся индексами узла. Оператор интерполяции с сетки p-1 на сетку p строится в виде произведения одномерных операторов интерполяции, перестановочных между собой:
Q =П ip] .
j=3
В этом случае значения сеточных функций на уровне разбиения p вычисляются по простым формулам, при этом на этапе интерполяции по 2-му и в особенности по 3-му направлениям выполняются одинаковые действия со всеми элементами больших массивов, которые реализуются с помощью специальны« программ на языке Assembler.
Аналогично выполняется операция сужения с сетки уровня p на сетку уров-
ня p-1, представимая в виде
qp=п .
j=1
Неточности, к которым в ряде случаев приводят интерполяция и сужение по целочисленным координатам, не страшны, посколыку в данном случае речь, идет о построении оператора близкой задачи.
Сужение матрицы Ар+1 на сетку уровня р дает матрицу Ар со сходной структурой. При этом Ару = 0"р+1 А^^^ Введем единичный вектор "Ь/р такой, что компоненты
(Я) = §(/1 Я )§(/2 Я )8(/з Яз ) , где Ъ(ш, п) - символ Кронекера; / = (/1, /2, /3) и я = (яь Я2, Яз) - мулытииндексы уз-
лов сетки уровня р.
Вектор Ару Ц является к-м столбцом блока Ару,
где
k = (/3 — 1)Mp + (¡2 — 1)MM2 + ¡1. Нетрудно показаты, что существует возможность. одновременного вычисления элементов многих столбцов блока Api с помощью 27 векторов
°p, Яъ Чз=12 3,
k
где k = (ki, k2, kp) и принимает те значения, для которых / пробегает по всем узлам сетки уровня p, при этом / = (q1 + 3(k1 — 1), q2 + 3(k2 — 1), q3 + 3(k3 — 1)).
Совокупность векторов Apij<5pq содержит все элементы ненулевых диагоналей Apij, которые помещаются на нужные места с помощью несложного алгорит-
ма после выполнения операции Qp+! Ар ^
были рассмотрены ранее.
Для формирования матрицы, соответствующей исходной аппроксимации задачи, используются 8-узловые конечные элементы, не требующие численного интегрирования. Конечный элемент строится путем осреднения по двум аппроксимациям тетраэдрами, как показано на рис. 2. В этом случае нетрудно получить формулы для вычисления вклада тетраэдров данной ориентации во все элементы ненулевой диагонали блока Аг;-, для чего также используются специаль-
б
а р
р +1CT q >
составные части которой
Рисунок 2
ные программы, написанные на языке Assembler.
Изложенная методика и алгоритмы были впервые реализованы авторами более 20 лет назад в виде программного комплекса РКТП (расчет конструкций в трехмерной постановке) по решению пространственных задач на ЭВМ класса ЕС [1]. Основные характеристики данного пионерного программного комплекса таковы. Большинство программных модулей было написано на языке PL/1. Для реализации арифметических операций с большими массивами чисел, для решения системы уравнений прямым методом, для быстрого обмена с внешней памятью (накопители на магнитных дисках) и для ряда других операций использовались программы, написанные на языке Assembler. Хорошая сходимость метода для широкого диапазона задач позволила сделать вариант программы промышленного типа, рассчитанный на широкого пользователя, от которого не требовалось квалификации вычислителя. Был разработан и внедрен ряд сервисных возможностей по вводу и контролю исходной информации. Использование вспомогательной крупной сетки позволило сделать удобным ввод координат узлов. Как правило, сначала вводились координаты узлов крупной сетки, затем, происходила генерация с помощью линейной интерполяции координат узлов исходной сетки, а затем, если это требовалось, осуществлялось исправление отдельных узлов или групп узлов исходной сетки.
Уже в первых версиях программного комплекса была реализована возможность автоматического использования всей доступной оперативной памяти ЭВМ, увеличение которой позволяло расположить большое количество массивов с размером, равным количеству узлов сетки, тем самым сокращая количество обменов и уменьшая время решения задачи. Была предусмотрена возможность поэтапного решения задачи. Каждый этап имел свою точку входа. В отдельные этапы выделялись ввод и проверка исходных данных, формирование исходной разрешающей матрицы, сужение исходной матрицы на все уровни сеточного разбиения, обращение матрицы на самой крупной сетке, итерационный процесс, дифференцирование решений, выдача результатов.
В последующих версиях программного комплекса реализовывался вариант, в котором матрицы Dp являлись диагоналями соответствующих матриц Ар (см. пункт 2). Были задействованы задачи теплопроводности, изотропной теории упругости и ортотропной теории упругости. Имелась возможность подключать другие виды краевых задач.
На рис. 3 демонстрируются сеточные разбиения реальных краевых задач, решение которых было осуществлено с помощью первых версий программного комплекса на ЭВМ ЕС-1033 в первой половине 80-х годов прошлого столетия. Все задачи имеют краевые условия, соответствующие их функциональному назначению. Критерий сходимости равен 8г- = (\\Аик -/\\Аи0 -/||г)-100%, I = 1,2,3, соответственно по нормам, рассмотренным в пункте 3. Все компоненты вектора начального приближения и0 равны 0, за исключением тех, которые соответствуют условиям 1-го рода.
На рис. 3 а показано сеточное разбиение при расчете теплового поля 1/16 симметричной части корпуса реактора. Решалась задача теплопроводности. Порядок системы превышал 30000. После 20 итераций 82= 1.71%, 83= 1.51%, 81 = 1.67%. Задача решалась по трехслойной схеме.
На рис. 3б показано сеточное разбиение при расчете фундаментной опоры. Решалась изотропная задача теории упругости. Порядок разрешающей системы уравнений превышал 7500. После 26 итераций 82 = 2.99%, 83 = 0.067%, 81 = 0.0054%. Задача решалась по двухслойной схеме.
На рис. 3в показано сеточное разбиение при расчете конструкции платформенного стыка крупнопанельных зданий. Решалась изотропная задача теории упругости. Порядок разрешающей системы уравнений превышал 33000. После 14 итераций 82 = 3.98%, 83 = 1.02%, 81 = 1.3 %.. Задача решалась по трехслойной схеме.
На рис. 3г показано сеточное разбиение при расчете фрагмента оболочки с цилиндрическим вырезом. Решалась ортотропная задача теории упругости. Порядок системы составлял 16000. После 8 итераций 82 = 0.93%, 83 = 0.37%, 81 = 0.86%. Задача решалась по трехслойной схеме.
Анализ результатов счета показал, что в подавляющем большинстве случаев при решении задач предлагаемым методом результат является достоверным при точности по \\ • \\2, равной 1 ^ 5%. При решении задачи, изображенной на рис. 3 а, по требованию заказчика результат посредством проведения дополнительных итераций был получен с точностью, превышающей 0.1%.
Литература
1. Белый М.В., Булгаков В.Е., Золотое А.Б. Полуитерационный многосеточный метод и его программная реализация для решения пространственных краевых задач. // Журнал вычислительной математики и математической физики, 1987, том 27, №6, с. 875-888.
2. Булгаков В.Е., Золотое А.Б., Белый М.В. Полуитерационный метод решения пространственных краевых задач расчета сооружений. // Строительная механика и расчет сооружений, 1985, №6, с. 38-40.
3. Золотое А.Б., Акимов П.А. Некоторые аналитико-численные методы решения краевых задач строительной механики: Монография - М.: Издательство АСВ, 2004. - 200 стр.
Рисунок 3
4. Золотое А.Б. К расчету трехмерных конструкций на ЭВМ. // Строительная механика и расчет сооружений, 1969, №6, с. 22-27.
5. Булгакое В.Е. К решению пространственной задачи теории упругости итерационным методом. // Численные методы и алгоритмы. М.: ЦНИИСК, 1981, с. 86-96.
6. Белый М.В., Булгакое В.Е. Алгоритмы автоматизированного решения пространственных краевых задач расчета сооружений полуитерационным методом. // Автоматизация расчета и проектирования промышленных и гражданских зданий и сооружений. М.: МИСИ, 1986, с. 112-123.
7. Федоренко Р. П. Релаксационный метод решения разностных эллиптических уравнений. // Журнал вычислительной математики и математической физики, 1961, том 1, №5, с. 1092-1096.
8. Бахеалое Н. С. О сходимости одного релаксационного метода при естественных ограничениях на эллиптический оператор. // Журнал вычислитель-
ной математики и математической физики, 1966, том 6, №5, с. 861-883.
9. Дьяконов Е.Г. Об использовании последовательности сеток при решении сильноэллиптических систем. // Вычислительные методы линейной алгебры. Новосибирск: ВЦ СО АН СССР, 1977, с. 146-162.
10. Марчук Г.И. Методы вычислительной математики. - М.: Наука, 1989. - 608 с.
11. Бахвалов Н.С., Жидкое Н.П., Кобельков Г.М. Численные методы. - М.: Лаборатория Базовых знаний, 2000. -624 с.
12. Федоренко Р.П. Введение в вычислительную физику. - М.: Издательство Московского физико-технического института, 1994. - 528 с.
13. Тимошенко С.П., Гудьер Дж. Теория упругости. - М.: Наука, 1975. -576 с.
14. Hachbusch W. On the convergence of miltigrid iterations. // Beitrage Numer. Math. 1981, Vol. 9, pp. 213-239.
