Разработка методов расчета задачи теплопроводности в контексте мультифизического подхода Текст научной статьи по специальности «Математика»

ЯлГЙОяО/

Уфа : УГАТУ. 2012____________________________^____________________________Т. 16, №3(48). С. 173-179

МАТЕМАТИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ ВЫЧИСЛИТЕЛЬНЫХ МАШИН

УДК 519.63

0. В. Шаповалов

РАЗРАБОТКА МЕТОДОВ РАСЧЕТА ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ В КОНТЕКСТЕ МУЛЬТИФИЗИЧЕСКОГО ПОДХОДА

В статье рассматривается вопрос разработки методов расчета теплопроводности в контексте совместного расчета многотель-ной модели и теплопроводности. Для использования в расчете задачи теплопроводности модифицирован и реализован метод Ньюмарка. Разработан метод для расчета жестких систем (в частности, задачи теплопроводности), обладающий высокой точностью расчета. Определена максимальная размерность задач, которые могут быть решены данным методом. Для расчета болыиихСЛАУ определенного вида был разработан алгоритм декомпозиции. Задача теплопроводности; междисциплинарный расчет; метод Ньюмарка; метод декомпозиции

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

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

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

Уравнения движения произвольной системы тел при таком подходе записываются в виде [5]:

Мх-ОтХ = /(ххЛ),

Вх = к(хл).

Здесь х - вектор переменных системы размерностью п, М - постоянная матрица коэффи-

Контактная информация: 8-919-980-56-09

циентов при вторых производных, Дх,х, 0 -вектор-функция правых частей дифференциальных уравнений, /) - матрица переменных коэффициентов уравнений связи размерностью к х п (к - число связей), И(х, х) - вектор правых частей уравнений связи, X - вектор множителей Лагранжа.

Уравнение является одной из форм записи уравнений Лагранжа 1 рода.

В ходе работы подвески жидкость в амортизаторе подвергается периодическому сжатию, что служит причиной повышения ее температуры. Совместное решение механической много-тельной модели и задачи теплопроводности носит итерационный характер. На каждом шаге расчета модели определяется положение амортизатора, на основе которого рассчитывается работа, совершенная демпфером. Затем, считая, что вся энергия идет на равномерный нагрев жидкости амортизатора, можно определить ее среднюю температуру, которая будет использоваться при решении задачи нестационарной теплопроводности в качестве граничных условий с внутренней стороны амортизатора. С внешней части граничным условием будет служить температура воздуха.

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

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

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

В процессе разбиения поверхность детали, представленная исходной САБ-геометрией, разделяется на полигоны. Разбиение поверхности на полигоны происходит с использованием свободной библиотеки функций для работы с геометрией ОрепСА8САББ. Также на первом этапе происходит разбиение модели машиностроительного изделия секущими плоскостями. Полигоны используются при построении сечения: ищутся линии пересечения полигона и секущей плоскости. Затем происходит разбиение сечения на равномерную ортогональную сетку путем проведения взаимно перпендикулярных лучей и поиска точек пересечения этих лучей.

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

Дискретная математическая модель соответствует уравнению теплопроводности в частных производных:

^2" = х, У,2,{),

дt

(1)

где и - температура, q - коэффициент теплопроводности материала, у - теплоемкость, р - плотность, с начальными условиями

и (х, у, z,0) = / (х, у, 2) (2)

и граничными условиями

и(х,у, 2,0 |г = ть. (3)

В качестве методов решения данной задачи использовались различные подходы. В качестве первого подхода рассматривалась простейшая разностная схема на основе явного метода интегрирования Рунге-Кутты четвертого порядка.

Уравнение (1) можно с заданной точностью представить в виде

д2

ди

аТТ + КРт- - ду и = в(х,У,2,*). (4)

д^ д*

При а ^ 0 решение уравнения (4) стремится к решению уравнения (1). Необходимое условие для величины а:

а<

у 2Р2

4Ч1'п

(5)

аТи(0 + УрГ(0 + дХ\Тп(*) = 0; п = 1,2,..., (6)

где Хп - собственное значение, Tn(t) - некоторая функция времени, соответствующая собственным значениям Хп. При достаточно малых а решение уравнения (6) стремится к решению уравнения

УрТп (*)+ Ч1\Тп ()= 0. (7)

Дискретизация уравнения (4) по методу конечных элементов приводит к уравнению вида

МХ + СХ + КХ = 0(х, t), (8)

где X - вектор узловых переменных размерности п, Мпхп - матрица при вторых производных, Спхп - матрица демпфирования, Кпхп - матрица жесткости, G(x, *) - вектор узловых нагрузок, соответствующий источникам тепла.

Матрица жесткости в предлагаемом методе складывается из элементов, которые для связи двух соседних узлов записываются в виде

к = д /1/9, (9)

где I - расстояние между соседними узлами.

Матрица демпфирования складывается из элементов вида

С =ур /1/9. (10)

В рассматриваемой реализации метода предполагается, что решетка разбиения - кубическая, с кубиками одинакового размера с длиной ребра I. Использование равномерной кубической решетки позволяет интегрировать уравнение (8) без явного вычисления матрицы жесткости. Правая часть уравнения (8) для -го узла будет вычисляться по формуле

т . . т . ,

ах, = -Ек(х< - хз)-Еси(х - хJ)+ 8(х*)’ (11)

3=1

В методе Фурье уравнению (4) соответствует уравнение для функции времени

где т - число контактирующих с узлом элементов.

Таким образом, был реализован модуль расчета задачи теплопроводности с использованием явной схемы интегрирования. Описанная выше технология расчетов позволяет использовать САБ-геометрию любой сложности, имеет достаточно высокую точность расчета, предоставляет высокую скорость разбиения геометрии на сетку, генерации и интегрирования дифференциальных уравнений. Выбранные методы разбиения модели и проведения расчетов допускают эффективное распараллеливание на многопроцессорных вычислительных системах. В качестве возможных путей распараллеливания были рассмотрены различные технологии.

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

Т аблица 1 Результаты распараллеливания задачи теплопроводности

Количество

время расчета, с

узлов разбиения, млн. GpuTesla с 1060 1 ядро corei7 920 2,667ГГц 4 ядра corei7 920 2,667ГГц

0,67 1,38 4,13 2,74

1,02 1,59 5,7 3,69

1,3 1,7 7 4,39

1,7 2 8,7 5,3

2,2 2,36 11,4 6,65

3,05 2,8 15 8,5

4,33 3,72 19,8 11,8

6,44 4,76 28,3 15,9

Визуализация результатов интегрирования дифференциальных уравнений происходит с использованием библиотек ОрепСАЗСАБЕ иУТК.

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

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

Используемые в пакете ФРУНД явные метода интегрирования - метод Эйлера и метод Рунге-Кутта четвертого порядка - хорошо работают, когда скорости происходящих в системе реакций соизмеримы. В случае, когда скорость одного процесса значительно превосходит скорость другого, явные схемы перестают работать, адекватно определяя поведение только одного из происходящих процессов. Системы обыкновенных дифференциальных уравнений, описывающие совокупности таких процессов, называются жесткими. Известно, что жесткие СОДУ эффективно решаются только с помощью неявных схем интегрирования.

Frame; 25

Рис. 1. Тепловой расчет гидропневматической рессоры

В качестве одного из неявных методов для реализации был выбран алгоритм Ньюмарка, который, правда, используется только в МКЭ-системах, но может быть использован после определенной модификации. Чтобы не отходить от терминологии метода будем обозначать изменение температуры перемещением, скорость изменения - скоростью, а ускорение изменения температуры - просто ускорением. Метод основан на предположении о линейном законе изменения ускорения w(t) на интервале |

/, + /| (рис. 2). Функции перемещения И'(/) и скорости w(t) произвольной узловой точки, совершающей колебания, представляются в виде отрезков степенных рядов

е

w{t,+t) =

w(ti+t)=wi+t-wi + a2t2wi

w.+t-wA------w.+a, tw.,

i i i i i ?

где a\,a2 - числовые параметры, определяющие остаточные члены степенных разложений функций и (/) и и’(/) на интервале | /, + /|.

[/„ /, + /]

Заменяя м’ приближенным разностным отношением (м/м - и’()//. преобразуем выражения к виду

+ а, ? (#,.+1 - #,),

мл+1 =м,1+1-й1 + а2 Найдем из первого равенства величину м>м :

^,+1 = л (71 Ґ

і і Г і л

-т(^і-------------------:Ч + (12)

(71 Ґ

2а

і /

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

™м=-{™м-™,) + а1 ґ

Ґ Ґ ( „ А

1- 2_ ч + +— 2- 2_

V аи 2 V

М/;.

(13)

Последние два выражения являются рекуррентными соотношениями метода Ньюмарка. Для вычисления векторов-столбцов перемещений {IV. скоростей {ж}; и ускорений

в момент tj+t = tм уравнение движения записываем в виде

М{ж}м+с{ж}м+к{ж}м ={Р(0},+1. (14)

После подстановки в него выражений (12) и (13) и преобразования к компактной форме, получим матричное уравнение вида

А{^}М={Р}М, (15)

где

А = К + ^С + —,

М+1=И0},+1 +

+м

Лм+^-М,+ т—іИ

а1 ґ а1 ґ 1 2а1 )

+

+ С

—м,+ --і И,

+ -

Сл

В начале алгоритма Ньюмарка, исходя из заданного шага интегрирования t, формируется матрица А, после чего она приводится к треугольному виду и до завершения вычислительного процесса не изменяется. В данном случае существует особенность

ф}м=шм. (16)

Сопоставляя формулы (15) и (16), получим

{5(01. = {5(01. - ф}+. - к{ж}м,

{5(0},+> = {5(01,+. + Ф1+. + к{ж}м.

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

Векторы-столбцы узловых ускорений и скоростей вычисляем по формулам:

М+.=4^}'+, -М)-+[г1-у-1и=

а1 Г а1 ? ^ 2 а1)

М,+1=М,+[(1 -сф\+«МЛ-

Для обеспечения устойчивости алгоритма прямого интегрирования параметры а\ и а2 должны удовлетворять следующим условиям: а2 > 0 и ах > 0,25 (а2 + 0,5 )2. Без учета искусственного рассеивания энергии максимальная точность метода Ньюмарка достигается при ах =

= 0,25 и а2 = 0,5.

Для запуска алгоритма метода Ньюмарка требуется задать в момент времени ( = 0 векторы-столбцы {Ж}0, [IV },. {ж }0 (начальные условия). В начальный момент времени будем считать скорости и ускорения нулевыми, а перемещения (температуры) должны задаваться. По умолчанию температура будет равна температуре окружающей среды.

Ключевым моментом решения является решение СЛАУ, после сравнения различных реализаций решение СЛАУ было реализовано с использованием библиотеки МКБ, что, с одной стороны, позволило сделать расчет быстрым и точным, но в то же время ограничило максимальный размер решаемых задач. Для реализации потребовалось вычислять матрицы жесткости и демпфирования и хранить в явном виде.

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

Для преодоления ограничения на размерность было принято решение разработать алгоритм декомпозиции [3]. Задача теплопроводности относится к задачам с одной степенью свободы, так как рассчитываемая система тел является жесткой, и структура не меняется со временем. Таким образом, было необходимо разработать новый метод расчета жестких систем. В общем случае, система тел описывается следующей системой Лагранжа 1 рода:

Мх - БТХ = /(х,х,Д £>х = к(х,х).

Для жесткой системы связи имеют более простой вид, и записывается следующим образом:

Мх- /)тХ = /(х,х,0,

Ох = 0.

Особенностью жесткой системы является ее

полносвязность, в связи с чем, решая систему все х оказываются одинако-

относительно

к

выми, т. е. V/ X. =//N, где х - вектор ускорений, а N — число узлов разбиения. Это знание можно использовать для более эффективного решения системы и решать ее не относительно

, а относительно X, т. е. подставляя найден-

ої ы

ные х в первое уравнение, получаем /У X = Мх - /(х. х. /). Данная система для абсолютного большинства случаев является вырожденной, т. е. не имеет решения классическими методами Гаусса или Крамера. Для таких систем существуют алгоритмы решения, наиболее подходящим для данного случая является алгоритм з\с1-разложсния.

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

В данной работе метод рассматривался применимо к тепловой задаче, х соответствует средней скорости нагревания тела, а X - тепловому потоку между соседними узлами.

У предлагаемого метода есть свои достоинства и недостатки. К достоинствам однозначно можно отнести большую точность и скорость расчета, поскольку метод является абсолютно устойчивым, так как погрешность не накапливается в связи с тем, что распределение температуры считается только в момент записи результатов. Этим же объясняется и высокая скорость расчета, поскольку время расчета не зависит от параметров расчета (времени интегрирования, шага расчета), а зависит от количества кадров времени, которые нужно записать. Но у метода также есть существенные недостатки: э \с1-

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

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

ОтХ =Мх - /(х,х^) = В.

Примем X = 1)у. тогда система примет следующий вид:

= В.

где /У/) - матрица жесткости системы. Теперь необходимо предобусловить матрицу, для этого достаточно добавить ненулевой элемент (от 1 до 100) к одному из диагональных элементов матрицы /У/). После этого система будет хорошо обусловленной и будет решаться стандартным методом Гаусса.

Эта методика была реализована, для решения СЛАУ использовалась функция РагсНво для решения СЛАУ для разреженных систем из библиотеки МаЙ1Кете11лЬгагу. Полученный метод обладает безусловной устойчивостью, высокой скоростью и точностью. Недостатком метода является ограничение на размерность решаемой системы, поскольку РагсПзо не решает системы с количеством ненулевых элементов больше 10 млн.

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

менее 1 %), при этом требуется порядка 8 Гб оперативной памяти. Увеличение размерности до 15 миллионов ненулевых элементов дает погрешность 7 %, что уже неприемлемо для большинства задач моделирования. Поскольку существует необходимость работы с большими массивами данных, надо каким-либо образом оптимизировать процесс решения СЛАУ с целью уменьшения требуемого объема памяти и времени решения. Одним из таких методов является замена одного решения СЛАУ сверхбольшой размерности [2] на решение нескольких СЛАУ меньшей размерности, основанное на знании структуры портрета матрицы.

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

Мх - И7X = /(ххл),

Вх = к{хл).

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

_ I1

------>

Рис. 3. Цепная система тел со связями

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

Л*; х2 х3 Р ♦< >♦< >♦ >

Я, Л-2

Рис. 4. Граф взаимодействия узлов в теле

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

хх + Х1 = О,

Х2 - А,; + Х2 = О, х3-Х2 = Р, х1-х2 = О, х2 - х3 = 0.

Та же система в матричной форме будет иметь следующий вид.

(1 0 0 1 0 " Г°1

0 1 0 -1 1 *1 0

0 0 1 0 -1 Хі = р

1 -1 0 0 0 к 0

.0 1 -1 0 ІА, л

Из графа видно, что каждый узел связан максимум с двумя соседними, а крайние узлы связаны только с одним телом. Если начать рассчитывать поведение модели с одного из крайних узлов, можно последовательно исключить все неизвестные системы. Таким образом, общее решение системы будет состоять из прямого хода, где на каждом шаге решение будет выражаться через одну переменную 7ч, и обратного хода. Обратный ход начинается после того, как будет решена последняя СЛАУ прямого хода, решение которого уже не будет зависеть ни от каких переменных. Во время обратного хода происходит последовательная постановка уже рассчитанных значений, данный этап полностью соответствует соответствующему этапу алгоритма Гаусса.

Для осуществления одного шага прямого хода необходимо решить систему из двух уравнений.

X, =°,

=0.

Поскольку зависимость ДЛЯ переменных Xj _ 1 и установлена на предыдущем шаге (для первого шага Х\ = Ю-\) = - Х}). подставляя выражение для х,- _ 1, получаем систему из двух уравнений и трех неизвестных, которое можно записать следующим образом:

(-1Л 1

[с і

где А - матрица подсистемы, С - константа, полученная сложением силы, приложенной к данному узлу, и постоянной части решения из предыдущего шага.

Пользуясь свойством линейности,

Ґ-О

С

= А-1 ■

- X, А-1 ■

Л

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

и решение для

( X ^

N

К-1

будет определено одно-

значно. Обратный ход представляет собой подстановку и достаточно тривиален.

Был проведен ряд тестов для цепочек разной длины, результаты тестов приведены в табл. 2.

Таблица 2

Расчет протяженного тела

Размерность Погрешность, % Время расчета, с

1,0Е+05 0 0,09

1,0Е+06 1,0Е-04 0,89

1,0Е+07 1,0Е-02 8,9

1,0Е+08 3,42 89,6

1,0Е+09 12,73 894

Для сравнения была выбрана самая быстрая на сегодняшний момент библиотека для решения СЛАУ для разреженных матриц 1п-1е1МКЬРагсПзо [4]. Для решения СЛАУ сверхбольшой размерности ей требуется значительный объем памяти, поэтому тесты проводились на вычислительном сервере кафедры ЭВМ и систем Волгоградского государственного технического университета, обладающем 16Гб оперативной памяти и процессором 1гЛе1Соге17 720. Данного объема памяти достаточно для решения системы размерностью пять миллионов.

Таблица 3 Расчет цепной системы в МКЬРапИво

Размерность Погрешность, % Время расчета, с

1,0Е+05 0 0,964

1,0Е+06 1,0Е-04 9,46

3,0Е+06 0,3 29,27

5,0Е+06 6,01 51,006

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

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

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

1. Антонов А. С. Параллельное программирование с использованием технологии ОрепМР: учеб. пособие. М.: Изд-во МГУ, 2009. 77 с.

2. Баландин М. Ю., Шурина Э. П. Методы решения СЛАУ большой размерности: учеб. пособие. Новосибирск: НГТУ, 2000. 76 с.

3. Банах Л. Я., Горобцов А. С., Чесно-ков О. К. Условия разбиения системы дифференциально-алгебраических уравнений на слабосвязанные подсистемы / Л. Я. Банах // Журнал вычислительной математики и математической физики. 2006. Т. 6, № 12. С. 2223-2227.

4. Библиотека РАМЖО. [Электронный ресурс]. Режим доступа: http://www.pardiso-project.org. 2011.

5. Виттенбург И. Динамика систем твердых тел. М.: Мир, 1980. 294 с.

6. Официальный сайт пакета ФРУНД [Электронный ресурс]. 2012. Режим доступа: Ьйр:/Лгшк1. vstu.ru.

7. Центр разработки ]\ГУГО1А С1ЮА. [Электронный ресурс]. Режим доступа: Ьир://^^\¥.пу1сИа. сот/оЬу ect/cuda_home_new. Мт1.

ОБ АВТОРЕ

Предложенный метод расчета имеет более высокую скорость расчета по сравнению со

Шаповалов Олег Владимирович, асп. каф. ЭВМ и систем Волгоградск. гос. техн. ун-та.