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

УДК 621.321

Разработка модели термодинамической системы на основе теории подобия и метода Монте-Карло в технологии параллельных вычислений CUDA

Д.М. Севрюгов, А. И. Тихонов ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,

г. Иваново, Российская Федерация Email: ait@dsn.ru

Авторское резюме

Состояние вопроса: В связи с ростом популярности технологии параллельных вычислений CUDA на многопроцессорных видеокартах фирмы nVIDIA значительные перспективы в плане решения задач теплообмена имеет метод Монте-Карло. Одним из наиболее естественных решений при моделировании термодинамической системы методом Монте-Карло выглядит подход, основанный на положениях молекулярной физики и теории подобия. Материалы и методы: Модель термодинамической системы основана на версии метода Монте-Карло и реализована путем имитации поведения молекул газа в соответствии с положениями молекулярной физики и теории подобия с применением технологии параллельных вычислений CUDA на графических процессорах nVIDIA. Результаты: Разработан математический аппарат имитационной модели термодинамической системы, адаптированной к технологии параллельных вычислений CUDA.

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

Ключевые слова: молекулярная физика, молекулярно-кинетическая теория, метод Монте-Карло, теория подобия, технология параллельных вычислений CUDA.

Development of a thermodynamic system model based on similarity theory and Monte Carlo simulation in CUDA parallel computing technology

D.M. Sevrugov, A.I. Tikhonov Ivanovo State Power Engineering University, Ivanovo, Russian Federation

Abstract

Background: As the CUDA technology of parallel computing on multi-processor graphics cards produced by the nVIDIA company is gaining popularity, the Monte Carlo method is becoming more promising in terms of solving heat exchange problems. One of the most natural approaches to simulation of thermodynamic systems using the Monte Carlo method is the one based on the provisions of molecular physics and the similarity theory.

Materials and methods: The thermodynamic system model is based on a version of the Monte Carlo method and implemented by simulating the behavior of gas molecules in accordance with the provisions of molecular physics and the similarity theory, CUDA parallel computing on GPU nVIDIA.

Results: A mathematical tool of simulation model of a thermodynamic system adapted to the CUDA parallel computing technology has been developed.

Conclusions: The developed thermodynamic system model can be used to simulate the research of thermodynamic processes and cycles in ideal and real gas. The model can also be combined with a finite-element model of thermal field in solids and used in transformer heat calculation.

Key words: molecular physics, molecular-kinetic theory, computing technology.

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

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

Monte Carlo method, similarity theory, CUDA parallel

метода конечных разностей. Поэтому в настоящее время для решения подобных задач растет популярность метода Монте-Карло. Этому способствует рост возможностей компьютерной техники, в частности, появление доступных для широкого круга пользователей многопроцессорных ускорителей на графических видеокартах. Особенно популярной является технология CUDA компании nVIDIA, с помощью которой разработчик может распараллелить расчетный процесс на 1000 и более ядер [1, 2].

Одним из наиболее естественных подходов при моделировании ТДС с использованием метода Монте-Карло выглядит подход, основан-

ный на положениях молекулярной физики и [3, 4] и теории подобия.

Следует отметить, что теория подобия создавалась, в первую очередь, в приложении к физическому моделированию явлений в технических системах (основные положения теории подобия изложены, например, в [5-7]). С развитием компьютерной техники физическое моделирование было практически вытеснено математическим моделированием. Однако с появлением технологии имитационного моделирования, позволяющей имитировать в реальном времени поведение структурных элементов моделируемых систем, наметилась тенденция к повторному использованию отработанного аппарата теории подобия на новой основе.

В качестве примера рассмотрим имитацию текучей охлаждающей среды (газа) множеством структурных элементов (молекул), движущихся в вакууме по законам классической механики [8, 9]. При моделировании ТДС каждой молекуле можно поставить в соответствие некий виртуальный объект (частицу). Множество таких частиц, двигающихся в реальном времени в виртуальном пространстве, представляет собой имитационную модель ТДС.

Так как количество молекул газа N в ТДС оказывается чрезмерно большим, то рационально воспользоваться методами теории подобия. При этом количество частиц в модели N можно задать по своему усмотрению, так что N' << N . При этом масштаб по количеству молекул является безразмерной величиной:

mN = N'/ N,

(1)

оказывается одним из независимых масштабов.

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

mI = /' / /,

(2)

где I - метрический размер оригинала; I' - соответствующий ему размер в модели. Так как для построения модели используется библиотека OpenGL [10], отображающая линейные размеры, заданные в дюймах, а размеры оригинала задаются в метрах, то масштаб по длине имеет размерность дюйм/м.

Время в модели удобно измерять не в секундах, а во фреймах, т. е. в кадрах изображения. При этом необходимо ввести коэффициент перевода, фрейм/с:

Kt = t ф / t c,

(3)

где t ф - время в модели, измеряемое во фреймах; t с - время в модели, измеряемое в секундах реального времени.

Масштаб по времени должен вычисляться с учетом коэффициента перевода:

mt = ^ = Kt — t t t t

где ^ - время оригинала, измеряемое в секундах.

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

Например, можно задаться желаемым значением наиболее вероятной скорости движения частиц и' в при какой-то характерной температуре Т (например, при Т = 273 К) и характерном давлении газа р (например, при р = 105 Па). Так, при и' в = 0,01 м/с движение частицы можно визуально проследить на экране монитора. При этом данной скорости частиц соответствует наиболее вероятная скорость молекул газа, м/с,

2kT

m,

2RT

0 V H

-23

(5)

где к = 1,3810- Дж/К - постоянная Больцмана; т0 - соответствует массе молекулы; Р = 8,31 Дж/(моль-К) - универсальная газовая постоянная; ц - молярная масса газа.

Скорость частиц в модели должна быть переведена из привычных единиц измерения м/с в единицы измерения модели дюйм/фрейм. При этом масштаб по скорости вычисляется как, (дюймс)/ (фреймм),

£в. m

ив Kt

(6)

(4)

Тогда масштаб по времени найдем как

т1 = т1 / ть. (7)

Поведение частиц модели, являющихся прообразами молекул газа, в общих чертах повторяет поведение молекул. Однако это подобие является неполным. В частности, отсутствует факт взаимного соударения частиц модели. Вместо этого каждая ¡-я частица совершает броуновское движение с длиной свободного пробега /' ¡, которая вычисляется с помощью генератора случайных чисел, возвращающего случайную величину 0 < ю ¡ < 1. При этом

/) =- 1п(ю,) •( / ■), (8)

где (/') - средняя длина свободного пробега частиц, выбираемая из условия соответствия ее средней длине пробега молекул газа (/). Масштаб по средней длине пробега молекул

тщ = Ц V <0 (9)

удобно выбрать в качестве четвертого независимого масштаба подобия, наряду с тд/, т/ и т

m =

Для этого необходимо задаться желаемым значением средней длины свободного пробега частиц модели (/') при какой-то характерной

температуре Т (например, при Т = 273 К) и характерном давлении газа р (например, при р = 105 Па), затем рассчитать среднюю длину пробега молекул газа при этой температуре по формуле 3п

</>=

пти

(10)

и по (9) вычислить т^. В (10) т0 - масса молекулы; п - динамическая вязкость газовой среды при заданной температуре Т и заданном давлении р; п - концентрация молекул газа при заданном давлении р: 3р .

п=

т0и

(11)

и - среднеквадратичная скорость молекул газа при заданной температуре Т:

и=

3кТ т0

(12)

к = 1,3810-23 Дж/К - постоянная Больцмана.

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

Скорость движения частиц модели на прямолинейном участке своей траектории выбирается исходя из температуры газа Т (в начале эксперимента температура газа принимается равной некоторому заданному начальному значению Т0). Скорость /-й частицы ина линейном участке траектории рассчитывается в соответствии с распределением Максвелла:

2

/ хЭ/2 три2

р(и)=ч 2тт] л~2кт. <13>

При этом реализуется следующий алгоритм. Для заданной температуре Т по (13) строится кривая распределения Максвелла (рис. 1).

Функция распределения скоростей

. 0.025

60 ВО V (скорость, м^с)

Рис. 1. Распределение молекул по скоростям при N = 30000 и Т = 100 К: 1 - полученное из модели для водорода; 2 - рассчитанное по (13)

По данной кривой строится обратная функция и(ю), представляющая собой зависимость

скорости частицы и от некоторой генерируемой случайным образом величины ю, значение которой изменяется в пределах от 0 до 1 (рис. 2).

Функция для определения скорости частицы

НО 130 120 110 100 ЭО ВО 70 60 50 40 30 20 10 О

р (случайное число)

Рис. 2. Функция для определения скорости молекул по результатам генерации случайного числа р

При вычислении скорости /-й частицы при заданной температуре Т сначала генерируется случайное число 0 <ю, < 1, а затем используется формула

и'/ = тии(ю,-). (14)

Полученное таким образом распределение скоростей частиц в расчетной области модели в состоянии термодинамического равновесия соответствует распределению Максвелла.

При столкновении частиц с изотермической стенкой происходит теплообмен. При этом скорость частицы пересчитывается по кривой Максвелла, соответствующей температуре стенки. Таким образом, изменяя температуру стенок цилиндра, дожидаясь установления термодинамического равновесия и вычисляя из модели распределение скоростей молекул, можно оценить, насколько велика расчетная погрешность, вызванная дискретностью модели. В частности, на рис. 1 приведены две кривые распределения молекул водорода по скоростям при Т = 100 К при количестве частиц в модели N = 30000. Реальное распределение отличается от вычисленного по (13). Однако с ростом количества частиц до нескольких миллионов эта погрешность снижается.

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

Я = т0

= т„

N

—— X Г(*

3^ 'АЭ' ^ Г\

] .пост+

) (в у.пост-)] :

/

6t' АЭ

N >

'X[т о (и'2 +-и'2-) ]=1

(15)

где Nt - количество частиц, ударившихся об элементарную поверхность стенки АЭ' за некоторое время t; / - количество степеней свободы моле-купы газа; (в'

].пост + /> ].пост-

)- энергия поступательного движения ]-й частицы до (индекс «-») и

после (индекс «+») столкновения со стенкой; и'/+ и и' у _ - скорость ./-й частицы до и после столкновения со стенкой.

Масса частицы т 0 соответствует массе молекулы т0 и определяется из условия равенства суммарной массы всех молекул ТДС суммарной массе всех частиц модели: т '0 N' = т0М , (16)

т. е.

, N

m 0 = — m0 = -

т. е.

N '

1

m

1

mN

-mo = mmmo,

(17)

(18)

где тт = 1/mN - масштаб по массе молекул.

Масштаб по удельному тепловому потоку тц определяется из выражения

q' mw q

qo =—,=z:?—ti=mqq ,

S ' t ' m;2mf St

т. е.

m,

mq =Za

w

mf 1 2

— m^ — mt

(19)

(20)

где масштаб по энергии тМ = т;2/ тг2 вычисляется исходя из выражения для внутренней энергии газа

М' = - vR 'Т' = тктт— vRT = 2 R Т 2 2 2 = М = т2 М = .

(21)

w

-n-Î mf

В свою очередь, масштаб по температуре

mT =

m

т =-2 вычисляется исходя из выражения

mmmt

для температуры газа, которое для ¡-го элементарного объема имеет вид

p -, = mp Pi = mf 1 pi

t ' = i = p ri = ' _ i, _ ~ „ /, ~

n 'i k mn n,k mN m;mf2 nk

(22)

= m; Pi =

mNmf2 nk

2 =

масштаб по давлению газа mp =

вычис-

m,m

;'"t

ляется из основного уравнения молекулярно-кинетической теории идеального газа, которое для ¡-го элементарного объема имеет вид

Р'/ = | Л '/ т 0 (и '2 ) = тпттти 1 п/т0 (=

mN 1 т2 1

= з ~--2 Р< =-1Р< = трР/,

т;3 mN т2 т/т2

масштаб по концентрации газа тп = т;у /т3 вычисляется исходя из выражения, которое для /-го элементарного объема имеет вид

(2f)

, = N 'i = -nN,- =

n , =

i ^ mf^ mf

n = mnni ,

(24)

где N, - количество частиц в /-м элементарном объеме модели V; N/■ - количество молекул в /-м элементарном объеме V.

Средний квадрат скорости частиц в /-м элементарном объеме модели в (24) вычисляется по формуле

1 N' .. N)

(и ?) = ^I и '2 = ^ До ^ +о % + и '2/), (25) / /=1 14 / /=1

где и 'х/, и 'у/, и '7/ - составляющие скорости /-й частицы в /-м объеме по осям координат.

Принимая ц' = ц, можно пересчитать универсальную газовую постоянную Р = 8,31 Дж/(моль-К), которая в модели принимает значение

п. k ц ' 1 k| 0 0

R = —î— =--— = mNR = mRR ,

m '0 mmm0

(26)

при этом mR = т;у .

Сила, действующая на поршень со стороны ТДС, вычисляется как сумма элементарных сил, действующих на каждую элементарную поверхность поршня площадью 5п:

F ' p 'i AS » mpm2 J pdSn

(27)

1

2c m, 2 m,F = —2F = mFF,

mm mf

где N'п - количество элементарных объемов, примыкающих к поршню в модели; S' /п - площадь контакта /-го элементарного объема модели с поршнем.

То есть масштаб по силе вычисляется как mF = mt / mf2 .

Работа, совершенная ТДС, равна механической энергии, переданной поршню:

Nx X2 m2

A • » £ F k Ax \ » mFm, J Fdx = -L A = mwA , (28)

k=1 x-, mt

где Nx - количество элементарных перемещений поршня; AX - элементарное перемещение поршня в модели.

Модель термодинамической системы была реализована в технологии CUDA на многопроцессорной видеокарте GeForce GTX 660 (960 ядер CUDA, тактовая частота 1032 МГцб, доступная графическая память 2 ГБ). В ТДС присутствовало 1000000 частиц. При этом скорость работы модели составляла 24 фрейма в секунду. Это позволяет достичь состояния термодинамического равновесия в ТДС при изменении температуры стенки цилиндра приблизительно за 10-30 с. Визуализация частиц в модели осуществляется с использованием библиотеки OpenGL [10].

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

п

1

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

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

1. Сандрес Дж., Кэндрот Э. Технология CUDA в примерах: введение в программирование графических процессов / пер. с англ. А.А. Слинкина. - М.: ДМК Пресс, 2011. -232 с.

2. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA. - М.: ДМК Пресс, 2010. - 232 с.

3. Гиббс Дж. Термодинамика. Статистическая механика. - М.: Наука, 1982. - 584 с.

4. Савельев И.В. Курс общей физики. Механика и молекулярная физика. - М.: Наука, 1982. - 511 с.

5. Веников В.А., Веников Г.В. Теория подобия и моделирования (применительно к задачам электроэнергетики): учебник для вузов по спец. «Кибернетика электрических систем». - 3-е изд., перераб. и доп. - М.: Высш. шк., 1984. - 439 с.

6. Веников В.А., Иванов-Смоленский А.В. Физическое моделирование электрических систем. - М.; Л.: Гос. энергетическое изд-во, 1956. - 359 с.

7. Иванов-Смоленский А.В. Электромагнитные поля и процессы в электрических машинах и их физическое моделирование. - М.: Энергия, 1969. - 304 с.

8. Севрюгов Д.М., Тихонов А.И. Разработка версии метода Монте-Карло для моделирования теплового поля // Состояние и перспективы развития электротехнологии (XVI Бенардосовские чтения): материалы Междунар. науч.-техн. конф. / Иван. гос. энерг. ун-т. - Иваново, 2011. -С. 134-137.

9. Свидетельство о государственной регистрации программы для ЭВМ «Виртуальный стенд-тренажер» № 2013611064 от 9.01.2013 / А.И. Тихонов, Д.М. Севрюгов. -М.: Федеральная служба по интеллектуальной собственности, патентам и товарным знакам, 2013.

10. Райт Р.С., Липчак Б. OpenGL. Суперкнига. - 3-е изд.: пер. с англ. - М.: Изд. дом «Вильямс», 2006. - 1040 с.

References

1. Sandres, Dzh., Kendrot, E. Tekhnologiya CUDA v primerakh: vvedenie v programmirovanie graficheskikh protsessov [CUDA technology in examples: introduction into GPU programming: translated from English by A.A. Slinkina]. Moscow, DMK Press, 2011. 232 p.

2. Boreskov, A.V., Kharlamov, A.A. Osnovy raboty s tekhnologiey CUDA [Fundamentals of CUDA programming]. Moscow, DMK Press, 2010. 232 p.

3. Gibbs, Dzh. Termodinamika. Statisticheskaya mek-hanika [Thermodynamics. Statistical mechanics]. Moscow, Nauka, 1982. 584 p.

4. Savel'ev, I.V. Kurs obshchey fiziki. Mekhanika i molekulyarnaya fizika [General physics study guide. Mechanics and molecular physics]. Moscow, Nauka, 1982. 511 p.

5. Venikov, V.A., Venikov, G.V. Teoriya podobiya i mode-lirovaniya (primenitel'no k zadacham elektroenergetiki) [Theory of similarity and modeling (in application to electrical engineering problems)]. Moscow, Vysshaya shkola, 1984. 439 p.

6. Venikov, V.A., Ivanov-Smolenskiy, A.V. Fizicheskoe modelirovanie elektricheskikh system [Physical modeling of electrical systems]. Moscow; Leningrad, Gosudarstvennoe energeticheskoe izdatel'stvo, 1956. 359 p.

7. Ivanov-Smolenskiy, A.V. Elektromagnitnye polya i protsessy v elektricheskikh mashinakh i ikh fizicheskoe modeli-rovanie [Electromagnetic fields and processes in electrical devices and their physical modeling]. Moscow, Energiya, 1969. 304 p.

8. Sevryugov, D.M., Tikhonov, A.I. Razrabotka versii metoda Monte-Karlo dlya modelirovaniya teplovogo polya [Development of a Monte Carlo method version for thermal field modeling]. Materialy Mezhdunarodnoy nauchno-tekhnicheskoy konferentsii «Sostoyanie i perspektivy razvitiya elektrotekhnologii» (XVI Benardosovskie chteniya) [Proceedings of the International scientific and technical conference "Status and development prospects (XVIth Benardos conference)]. Ivanovo, 2011, pp. 134-137.

9. Tikhonov, A.I., Sevryugov, D.M. Svidetel'stvo o go-sudarstvennoy registratsii programmy dlya EVM «Virtual'nyy stend-trenazher» № 2013611064 [State registration certificate of computer program «Virtual simulator stand» no. 2013611064]. Moscow, Federal'naya sluzhba po intellektual'noy sobstvennosti, patentam i tovarnym znakam, 2013.

10. Rayt, R.S., Lipchak, B. OpenGL. Superkniga [OpenGL. Superbook]. Moscow, Izdatel'skiy dom «Vil'yams», 2006. 1040 p.

Севрюгов Денис Михайлович,

ФГБОУВПО «Ивановский государственный энергетический университет имени В. И. Ленина», соискатель,

e-mail: zeron88@yandex.ru Тихонов Андрей Ильич,

ФГБОУВПО «Ивановский государственный энергетический университет имени В. И. Ленина», доктор технических наук, доцент, зав. кафедрой физики, e-mail: ait@dsn.ru