О приеме диагонализации матриц в методе конечных элементов при решении нестационарных задач

Размещено на http://www.allbest.ru/

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

О приеме диагонализации матриц в методе конечных элементов при решении нестационарных задач

Л.Н. Ясницкий

Аннотация

Излагается точка зрения автора на некоторые приоритетные вопросы в области развития и применения метода конечных элементов. Рассмотрен алгоритм применения метода конечных элементов с диагонализированными матрицами разрешающих СЛАУ в задачах термогравитационной конвекции жидкости. Приводится пример применения алгоритма для решения технической проблемы получения высококачественных стальных отливок методом математического моделирования.

Ключевые слова: метод конечных элементов; краевая задача; механика жидкости.

The point of view of the author on some priority questions in the field of development and application of final element method is stated. The algorithm of application of final element method with diagonalized matrixes of resolving systems of the linear algebraic equations in problems of heat convection of liquid is stated. It is given an example algorithm applications for the solution of a technical problem of receiving high-quality steel body by a method of mathematical modeling.

Key words: finite element method; boundary problem; mechanics of liquid.

Содержание

Введение

1. Алгоритм МКЭ с диагонализацией матриц СЛАУ в задачах термогравитационной конвекции жидкости

2. Пример математического моделирования

Заключение

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

Введение

Метод конечных элементов (МКЭ), как известно [1], был впервые предложен и применен в начале 40-х гг. XX столетия инженерами-механиками, представлявшими упругие тела в виде наборов брусьев и балок [2], и уже потом было получено его математическое обоснование. Во второй половине XX в. выходят основополагающие монографии О. Зенкевича, Р. Галлагера, Дж. Одена, Л. Сегерлинда, Г. Стренга и Дж. Фикса, Л.А. Розина и др. Методом конечных элементов решаются краевые задачи механики деформируемого твердого тела (МДТТ) в стационарной постановке. Применение же МКЭ для нестационарных задач МДТТ несколько задерживается. Причина в том, что в таких задачах требуется многократное решение разрешающих систем линейных алгебраических уравнений (СЛАУ), что на ЭВМ с ограниченным быстродействием требует большого расхода машинного времени.

Как выход из положения вводится упрощающая гипотеза - вместо центра тяжести масса каждого элемента равномерно распределяется по его узлам [3]. В результате матрицы разрешающих СЛАУ становятся диагональными. Необходимость решения СЛАУ отпадает и применение МКЭ в нестационарных задачах механики деформируемого твердого тела становится реальностью. Появляются публикации, в которых МКЭ успешно применяется в инженерных расчетах динамических систем и процессов (например, [4] и др.)

Надо отметить, что введение гипотезы равномерного распределения массы элемента по его узлам соответствует творческому стилю инженеров, предложивших и реализовавших вначале 1940 гг. саму идею метода конечных элементов [2]. Эта гипотеза понятна на интуитивном уровне, но с точки зрения классической механики массу элемента можно равномерно рассредоточивать по его узлам только в том случае, если элемент имеет правильную форму. В общем же случае для элементов неправильной формы (какие применяются в действительности) эта операция неправомерна и ее применение означает нарушение законов классической механики. Понятно, что такое нарушение вносит погрешность в результат решения краевой задачи. Но понятно и другое. Эта дополнительная погрешность, вызванная нарушением законов механики, стремится к нулю с измельчением конечноэлементной сетки. И значит, введение гипотезы равномерного распределения массы элемента по его узлам в методе конечных элементов вполне оправдано.

Тем временем вопрос о применении МКЭ для решения нестационарных краевых задач механики жидкости все еще остается нерешенным. В 70-х гг. делаются попытки применения метода конечных элементов для решения краевых задач механики жидкости в стационарной [5-7] и в нестационарной [8, 9] постановках, однако последние решаются только для существенно упрощенных условий, причем с нереально большими затратами машинного времени. Причина все та же. Матрица разрешающей СЛАУ, которая в нестационарных задачах должна решаться множество раз, имеет довольно широкую ленту, а значит, на каждом шаге по времени требуется применение ресурсоемких методов решения СЛАУ.

Автором настоящей работы в 1980 г. депонируется статья [10], в которой предлагается алгоритм диагонализации матрицы СЛАУ метода конечных элементов при решении нестационарных задач механики жидкости. Согласно этому алгоритму предлагается заменить каждый диагональный коэффициент суммой всех коэффициентов строки, которой он принадлежит, а коэффициенты, не лежащие на главной диагонали, - нулями. В результате отпадает необходимость решения СЛАУ, а значит, и применения ресурсоемких методов решения СЛАУ на каждом временном шаге. Машинное время, необходимое для решения краевой задачи, сокращается на несколько порядков.

Заметим, что с точки зрения чистой математики применение предложенного в работе [10] приема диагонализации матриц СЛАУ, на первый взгляд, кажется безграмотной затеей. Действительно, решения СЛАУ, полученные с помощью приема диагонализации матриц, в принципе не должны совпадать с решениями СЛАУ, полученными классическими способами, например, методом Гаусса, Крамера, Зейделя и др. Поэтому в упомянутой статье [10] в качестве обоснования правомерности применения приема диагонализации матриц приводится следующий довод. Во-первых, обращается внимание на то, что диагонализация матриц СЛАУ в задачах механики деформируемого твердого тела имеет физическую интерпретацию, а именно применение этого приема равносильно введению гипотезы равномерного рассредоточения массы элементов по его узлам. Как отмечалось выше, эта гипотеза не вызывает возражений, поскольку погрешность, обусловленная ее применением, стремится к нулю с измельчением конечноэлементной сетки.

В задачах же механики жидкости подобной физической интерпретации, к сожалению, придумать не удается. Поэтому в статье [10] высказывается предположение, что между задачами МДТТ и задачами механики жидкости существует некая математическая аналогия, а именно в задачах механики жидкости, как и в задачах МДТТ, погрешность, вносимая применением приема диагонализации матриц СЛАУ, с измельчением конечноэлементной сетки тоже стремится к нулю. И это предположение подтверждается путем проведения вычислительных экспериментов 10].

Правомерность и эффективность предложенного в работе [10] приема диагонализации матриц разрешающих СЛАУ подтверждается последующими публикациями автора настоящей статьи [11-27], посвященными применению метода конечных элементов для решения инженерно-технических задач, а также тем, что этот прием успешно запрограммирован в ряде современных пакетов прикладных программ (например ANSYS, LS-DYNA, FLUENT и др.), которые широко применяются во всем мире для решения различного рода задач, связанных с моделированием нестационарных гидродинамических явлений и процессов.

1. Алгоритм МКЭ с диагонализацией матриц СЛАУ в задачах термогравитационной конвекции жидкости

Покажем применение алгоритма диагонализации матриц разрешающей СЛАУ на примере применения метода конечных элементов для моделирования процесса получения стальных отливок.

Система уравнений термогравитационной конвекции вязкой несжимаемой жидкости в приближении Буссинеска [28] с введением в качестве единиц измерения длины - характерного линейного размера области , температуры - характерной разности температур , времени - , давления - записывается в безразмерной форме

(1)

(2)

(3)

Здесь введены числа Грасгофа

и Прандтля

,

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

Применяя схему искусственной сжимаемости [29] вместо (3) будем использовать уравнение

(4)

где - параметр релаксации.

Применение процедуры Галеркина к уравнениям (1), (2) и (4) с последующим преобразованием интегралов от вторых производных по теореме Грина приводит к уравнениям

(5)

(6)

, (7)

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

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

(8)

(9)

(10)

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

Выбирая в качестве весовых функций

и подставляя разложения (8), (9), (10) в (5), (6), (7), получаем

(11)

(12)

(13)

где

(14)

Здесь - номера узлов сеточной области.

Равенства (11)-(13) представляют собой систему обыкновенных дифференциальных уравнений первого порядка, в которых неизвестные функции зависят от времени. Если в этих уравнениях производные по времени представить в конечно-разностном виде, то интегрирование системы уравнений (11)- (13) сведется к многократному решению систем линейных алгебраических уравнений. Как уже отмечалось ранее, чтобы избежать этой ресурсоемкой процедуры, в работе [10] было предложено подвергнуть матрицу диагонализации путем замены каждого диагонального коэффициента суммой всех коэффициентов строки, которой он принадлежит, а коэффициентов, не лежащих на главной диагонали - нулями.

В результате такой процедуры уравнения (11)-(13) преобразуются к виду

(15)

(16)

(17)

где - символ Кронекера, - единичная матрица.

Последовательность интегрирования уравнений (15)-(17) принимается следующей. На -м шаге по времени из уравнения (15) находятся пробные скорости :

(18)

а из уравнения (16) - значения температуры:

(19)

Далее из уравнения (17) ищется поле давления, соответствующее пробным скоростям:

(20)

Здесь

,

- величина шага по времени, и - последовательные приближения к и . Процесс (18), (20) повторяется до тех пор, пока

.

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

Во всех решаемых по описанной методике задачах аппроксимации (8)-(10) заранее удовлетворяли граничным условиям задачи, поэтому в уравнениях (18)-(20) индекс граничные узлы области не охватывал. Если же для давления граничных условий в формулировке задачи не ставилось, то в уравнении (20) индекс охватывал и граничные узлы расчетной области и, таким образом, значение давления на границе получалось как естественное граничное условие.

2. Пример математического моделирования

На рисунке приведены результаты моделирования процесса затвердевания стальных отливок. Результаты получены путем совместного решения дифференциальных уравнений термогравитационной конвекции жидкости (1)-(3) и условия сопряжения теплового и гидродинамического полей при , в котором - температура выливаемости, соответствующая состоянию затвердевающего сплава, когда в единице объема содержится 30% твердой фазы.

Процесс затвердевания стальных отливок в моменты времени: а - 600 с; б - 900 с; в - 1800 с; д - 5400 с

Интегрирование уравнений (1)-(3) выполнялось методом конечных элементов с помощью алгоритма (18)-(20), реализующего идею диагонализации матриц СЛАУ [10]. Результаты представлены на рисунке в разные моменты времени: а, б,…, д в виде распределения изолиний функции тока в незатвердевших частях отливок (). Затвердевшие области отливок () на рисунке заштрихованы.

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

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

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

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

Заключение

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

Применение приема диагонализации матриц разрешающей СЛАУ метода конечных элементов для решения краевых задач механики жидкости позволило выполнить математическое моделирование ряда инженерно-технических процессов [11-27]. Этот прием, по-видимому, впервые предложенный автором настоящей статьи в работе [10], нашел применение в современных пакетах прикладных программ, таких как ANSYS, LS-DYNA, FLUENT и др. В настоящей статье изложен алгоритм его реализации на примере решения задачи выбора рациональных режимов получения стальных отливок методом математического моделирования.

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

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

2. Hrennikoff F. Solution of problems in elasticity by the framework method // Journal of applied mechanics. 1941. №8. P.169-175.

3. Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975. 320с.

4. Римм Э.Р., Нихамкин М.Ш. Об одном варианте метода конечных элементов для решения динамических задач теории упругости // Полимерные материалы в машиностроении: межвуз. сб. науч. труд. ПГУ и ППИ. Пермь, 1977. № 214. С.129-134.

5. Smith S.L., Brebbia C.F. Finite element solutions of Navier-Stokes equations for transient two-dimensional incompressible flow // Journal of computational physics. 1975. № 17. P.235-245.

6. Taborrok B., Lin R.C. Finite element analysis of free convection flows // Internation journal of heat and mass transfer. 1977. Vol.20, № 9. P. 945-952.

7. Oden J.T., Wellford L.C. Analysis of flow of viscous fluids by the finite element method // AIAA Journal. 1972. Vol.10, №12. P.1590-1599.

8. Kawahara M., Yoshimura N., Nakagawa K., Oshaka H. Steady and unsteady finite element analysis of incompressible viscous fluid // International Journal of Numerical Methods in Engineering. 1976. № 10. P. 437-456.

9. Usuki S. The application of variational finite element method to problems in fluid dynamics // International Journal of Numerical Methods in Engineering. 1977. №11. P.563-577.

10. Самойлович Ю.А., Ясницкий Л.Н. Алгоритм решения задач термогравитационной конвекции вязкой несжимаемой жидкости методом конечных элементов / Перм. гос. ун-т. М., 1980. Деп. в ВИНИТИ 24.03.80. №1131-80.

11. Самойлович Ю.А., Ясницкий Л.Н. Сопряженная задача теплообмена, гидродинамики и затвердевания расплава // Инженерно-физический журн. 1981. Т.XLI, №6. С.1109-1118.

12. Самойлович Ю.А., Ясницкий Л.Н. Математическое моделирование затвердевания стали с учетом термоконвективного движения расплава // Изв. вузов. Черная металлургия. 1981. №12. С. 75-78.

13. Самойлович Ю.А., Кабаков З.К., Ясницкий Л.Н. Сопряженная задача теплообмена и гидродинамики в затвердевающем расплаве // Теплофизика высоких температур. 1981. Т.XIX, №4. С.814-820.

14. Самойлович Ю.А., Емельяненко Ю.Г., Ясницкий Л.Н. Анализ завтердевания слитка ЭШП с учетом гидродинамических явлений в металлической ванне // Докл. АН УССР. Сер. "А". Физико-математические науки. 1981. №8. С. 91-94.

15. Самойлович Ю.А., Кабаков З.К., Ясницкий Л.Н. Математическое моделирование тепловых и гидродинамических явлений процесса затвердевания непрерывного слитка // Изв. АН СССР. Металлы. 1982. №2. С.62-68.

16. Медовар Б.И., Самойлович Ю.А., Емельяненко Ю.Г., Ясницкий Л.Н. Анализ тепловых и гидродинамических явлений в металлической ванне при электрошлаковом переплаве // Проблемы специальной электрометаллургии. Киев: Наукова думка, 1982. Вып.16. С.29-35.

17. Самойлович Ю.А., Ясницкий Л.Н. Неустойчивость тепловой гравитационной конвекции в жидком ядре затвердевающей отливки // Теплофизика высоких температур. 1982. Т.XX, №5. С.1002-1004.

18. Самойлович Ю.А., Кабаков З.К., Ясницкий Л.Н. Инженерная методика расчета параметров электромагнитного перемешивания расплава в условиях непрерывного литья // Изв. вузов. Черная металлургия. 1983. №10. С. 155-156.

19. Самойлович Ю.А., Кабаков З.К., Ясницкий Л.Н. Исследование термогравитационной конвекции при затвердевании жидкой стали методом математического моделирования // Инженерно-физический журн. 1983. №3. С. 456-473.

20. Самойлович Ю.А., Кабаков З.К., Ясницкий Л.Н. Гидродинамические явления при затвердевании непрерывного слитка в условиях индуктивного МГД-воздействия // Магнитная гидродинамика. 1983. №4. С.123-130.

21. Самойлович Ю.А., Кабаков З.К., Ясницкий Л.Н. Инженерная методика расчета электромагнитных перемешивающих устройств на машинах непрерывного литья // Магнитная гидродинамика. 1984. №2. С.120-126.

22. Федотов В.М., Субоч В.Д., Тихонов Н.И., Самойлович Ю.А., Ясницкий Л.Н. Физическое моделирование процесса кондукционного перемешивания расплава в непрерывном затвердевающем слитке // Магнитная гидродинамика. 1984. №4. С.95-100.

23. Самойлович Ю.А., Фультахт В.В., Кабаков З.К., Ясницкий Л.Н. Тепловая конвекция жидкой стали при затвердевании непрерывного слитка // Изв. АН СССР. Металлы. 1985. №1. С. 49-54.

24. Кирко И.М., Самойлович Ю.А., Долгих В.М., Хрипченко С.Ю., Ясницкий Л.Н. Электровихревой способ перемешивания расплава затвердевающих слитков // Магнитная гидродинамика. 1985. №3. С.100-107.

25. Медовар Б.И., Самойлович Ю.А., Емельяненко Ю.Г., Девингталь Ю.В., Андриенко С.Ю., Чайковский А.И., Чудновский А.Ю., Ясницкий Л.Н. Исследование электрошлакового процесса методами физического и математического моделирования // Проблемы специальной электрометаллургии. 1985. №3. С.5-10.

26. Андриенко С.Ю., Емельяненко Ю.Г., Ясницкий Л.Н. Исследование гидродинамики при электрошлаковом переплаве методом физического и математического моделирования // Проблемы специальной электрометаллургии. 1987. №4. С.5-7.

27. Андриенко С.Ю., Ясницкий Л.Н. Механизм удаления неметаллических включений при электрошлаковых процессах // Расплавы. 1990. №4. С. 114-115.

28. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972. 392 с.

29. Chorin A.J. Numerical solution of the Navier-Stokes equations // Mathematics of Computation. 1968. Vol.22, № 104. P.745-762.

30. Беляев Ю.П., Сердюков Г.В. Роль конвективных потоков в процессе кристаллизации крупных листовых слитков спокойной стали // Теория и практика металлургии чугуна и стали: Тр. Ждановского металлургического ин-та. Жданов, 1971. Вып.14. С.221.

31. Ефимов В.А. Теплофизические процессы и методы управления формированием отливок и слитков // Тепловые процессы в отливках и формах: Тр. XV совещания по теории литейных процессов. М.: Металлургия, 1972. С.15.

Размещено на Allbest.ru