Компьютерное моделирование процессов с распределенными параметрами Текст научной статьи по специальности «Математика»

УДК 519.8(004)

А. Л. Королев

КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ С РАСПРЕДЕЛЕННЫМИ ПАРАМЕТРАМИ

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

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

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

Пусть имеется теплообменник, представляющий собой канал, нагрев которого производится, например, путем пропускания электрического тока через его стенки (рис. 1). Примем следующие допущения: в теплообменник поступает холодная жидкость с температурой Т0; скорость движения жидкости внутри канала считается неизменной; плотность и теплоемкость теплоносителя слабо зависят от температуры и могут считаться неизменными. Температура теплоносителя по сечению канала считается одинаковой.

Ч

х

Рис. 1 Расчетная схема объекта исследования

В начальный момент времени ґ = 0 включается обогрев канала, т.е. q > 0 при ґ > 0. Здесь q - количество тепла, которое передается через стенки канала к теплоносителю в единицу времени через единицу площади поверхности канала.

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

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

( АТ /-)Т7 Л

51 рс— + сри — | = qП, 5 = кЯ2, П = 2кЯ,

^ дt дх)

где 5 - площадь поперечного сечения канала теплообменника; П - его периметр; р - плотность теплоносителя; с - его теплоемкость; и - скорость движения теплоносителя.

Физический смысл уравнения состоит в том, что его левая часть отражает изменение запаса тепла в элементе канала длиной Ах (рис. 2), а правая часть уравнения учитывает количество тепла, которое подведено к элементу канала через боковую поверхность за счет обогрева.

Рис. 2 Элемент канала или стержня

Простейшее преобразование исходного уравнения дает следующее:

дТ дТ qП „

— + и — = —— = О. дt дх Б • рс

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

- х - / Ь Л

виде: х = — I = t/1 — |.

Ь ^ и )

Преобразование уравнения дает соотношения

дТ и.дк=о- Ьдк дТ=Ок- дТ дТ=кО

дt Ь дх ; и дt дх и ’ дt дх и

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

зависимостью у = (Т - То)/1 —О|. Тогда безразмерное уравнение переноса,

Iи )

краевые и начальные условия примут вид

ду + ■ду = 1, у(^, х = 0) = 0, у(Г = 0,х) = 0.

дt д х

Обращаем внимание, что безразмерная модель не имеет параметров. Таким образом, все процессы переноса подобны между собой. Подробно метод построения безразмерных моделей описан в работе [5].

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

для дискретных моментов времени Xj±l = Xj ± h, ^+\ =^+г, h - шаг по ко-

ординате, т - шаг по времени.

Введем индексную форму обозначения переменных: у” = у(Уи, х^),

у”!+1 = У^„+\, Xj) , уП +! = У^п+\, Xj_1) . Пусть в момент времени t = tn известны значения уП в любой точке х}- канала. Конечноразностные соотношения для производных имеют вид

Последнее соотношение есть окончательная расчетная формула, применяя которую для j = 2 _ М (М - номер последней точки по координате х),

при х = 0, значение у ” всегда известно, т.к. задано краевыми условиями.

Результат решения задачи переноса средствами электронных таблиц в координатах у( х) для различных моментов времени представлен на рис. 3, 4. Электронная таблица по рис. 3 в целом содержит значения У0-У50. В ней

Моделирование процесса теплопроводности рассмотрим на примере следующей задачи. Имеется длинный металлический стержень (рис. 5), который в начальный момент времени имеет одинаковую температуру по всей длине Т0. В начальный момент времени температура на правом торце стержня (х = Ь) изменяется скачкообразно до значения Ті и сохраняет это значение. Левый торец стержня (х = 0) сохраняет постоянную температуру Т0 в любой момент времени. По боковой поверхности стержень взаимодействует с окружающей средой (нагревается или охлаждается).

Эх И дґ т

И позволяют представить исходное уравнение в следующем виде:

т

И

Преобразуем конечно-разностный аналог уравнения:

1

по известным значениям у ^ получим искомые величины Уj

п

Рис. 3 Электронная таблица решения задачи переноса

Рис. 4 Результат решения задачи переноса в среде электронных таблиц

Целью моделирования является определение закона изменения температуры стержня после скачкообразного изменения температуры его правого торца. Примем, что длина стержня значительно больше его радиуса Ь»Я. Материал стержня имеет достаточно высокую теплопроводность, следовательно изменением температуры в поперечном сечении можно пренебречь (она считается однородной в любом поперечном сечении). Таким образом, процесс распространения тепла по стержню можно описать одномерной моделью.

На рис. 5 представлена расчетная схема объекта исследования. Процессы передачи тепла в стержне описываются следующим уравнением:

ЭТ .Э 2Т + п Т Т )

р-с эТ+7 “(Т° Ь

где х - продольная координата; ґ - время; П - периметр поперечного сечения стержня; 7 - площадь поперечного сечения; а - коэффициент теплоотдачи в окружающую среду через боковую поверхность стержня; X - коэффициент теплопроводности материала стержня; р - его плотность; Т0 - температура окружающей среды; Т(ґ, х) - температура стержня - функция времени и координаты.

Ц >> Я

х = 0 Тп Ах х = Ц

Рис. 5 Расчетная схема моделирования процесса теплопроводности

В начальный момент времени Т(ґ = 0, х) = Т0, а краевые условия имеют вид

Т(ґ, х = 0) = То, Т(ґ, х = Ц) = Ть

Уравнение теплопроводности построено на основе закона сохранения энергии и закона Фурье. Закон Фурье связывает величину осевого потока тепла по стержню за счет теплопроводности с градиентом температуры:

q = -X—. Тепловой поток через боковую поверхность стержня определяется Эх

законом Ньютона: q = а(Т) - Т) . Данные законы имеют чисто феноменологический характер. В элементе стержня длиной Ах (рис. 2) протекают следующие процессы: накопление тепла, продольная передача тепла по стержню теплопроводностью и теплоотдача через боковую поверхность в окружающую среду.

Предварительно преобразуем исходную задачу:

ЭТ Э2Т , , > X , Па

— = а—— + к (То - Т), а =-----, к = .

дґ Эх2 10 ' р-с’ 7-р-с

В качестве масштаба по продольной координате выберем длину стержня Ц. Тогда безразмерная координата будет определяться соотношением х = х / Ц . После преобразования исходное уравнение примет вид

Ц ЭТ = Э2Т + кЦ Т - ) а ЭГ Эх2 а Т 0 ).

В качестве масштаба по времени примем t = Ь / а . Безразмерное вре-

— *

мя определим следующим образом: t = t/1 . Безразмерную температуру зададим соотношением у = (Т - То)/(Т - То). Уравнение теплопроводности, начальные и краевые условия в итоге преобразуются к следующему безразмерному виду:

ду =

ЭГ дх2

- ку, у ( = 0, х ) = 0, у (, х = 0 ) = 0, у (, х = 1) = 1.

- кЬ2

Безразмерный параметр к =------ является критерием подобия и харак-

а

теризует соотношение интенсивности процессов теплообмена с окружающей средой и теплопроводности по стержню. Заметим, что для значения к = 0 (теплоизоляция боковой поверхности стержня) в уравнении нет параметров, следовательно в этом случае для любых параметров все реальные процессы теплопроводности будут подобны.

Для численного решения задачи моделирования теплопроводности применим метод конечных разностей [2, 3]. Таким образом, решение будет найдено для конечного числа точек Ху в фиксированные моменты времени tn . Производные приближенно зададим соотношениями

у!+1 - у'\

ду

ді

п+1

п+1

п+1

2

. у]+1 - 2у] + у ]-1 _ д2у

дх 2

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

уП+1 уп уП+1 2 уп+1 + уп+1

у] - уу_уу+1 - 2 у] + у] -1 7-,п+1

- ку ]

1 п+1 | 1 + 2 .7 | п+1 + 1 п+1

ТУуп-1 -|Т + 7Г+к Іуп + 7Гуп+1

т

Последнее соотношение справедливо для значений у = 2 + (Ы -1)

(внутренние узловые точки по х). Для краевых точек значения упу+1 заданы

краевыми условиями: у = 1, уП+1 = 0 , а для у = М, уП+1 = 1. Значения уП в любой точке ху стержня в момент времени t = ^ считаются известными. Индексная форма записи переменных означает следующее: у’у = у^п,Ху),

уП+1 = y(tn+l, Ху), уп+1 = y(tn+l,Ху±1) и т.д., причем Ху±1 = Ху ± к, ^+\ = ^+ т.

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

уп+1 . Подобная специфическая система линейных уравнений решается с по-

т

мощью метода прогонки [2, 3], она может быть представлена в следующем каноническом виде:

Решение ищется в виде рекуррентного соотношения: у у = аууу+1 + в у .

Коэффициенты ау и в называются прогоночными. Положим у = 1, тогда у = а1 у2 + Р1, в то же время у = х у2 + ^1, из чего следует, что а = Х1, Р1 = ^1. Расчет остальных коэффициентов выполняется по формулам

После определения коэффициентов а у и ву решение выполняется по формуле у у = аууу+1 + ву , где у = (И -1) ь 2 .

Таким образом, алгоритм решения задачи включает: расчет значений коэффициентов Ау, В у, Су, ¥у и %1, ц>1, Х2 , ^2; задание значений коэффициентов а1, в1; расчет а у, в у , причем у изменяется от у = 2 до у = (М - 1);

расчет уГу+1, в этом случае у изменяется от у = (М - 1) до у = 2.

Результат численного решения задачи теплопроводности для различных моментов времени и определенного значения параметра к представлен на рис. 6, 7.

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

ду 0

разована для решения задачи с краевым условием второго рода: — = 0.

дхх=0

Это соответствует теплоизоляции левого торца стержня.

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

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

Аууу-1- Сууу + Вууу+1 =- ,

у1 =Х1 у2 + иЬ

уИ = х2 уИ -1 + и2.

В данном случае коэффициенты системы имеют значения

а у =----------------

у Су - А а у-1

Су - Ау а у-1

Рис. 6 Электронная таблица решения задачи теплопроводности

Рис. 7 Результат решения задачи теплопроводности

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

смотренных выше. Модель совместного протекания процессов строится также на основе закона сохранения энергии:

I dT + dT У ,dzT + п ! T ^ рсI — + u— I = Л—— +—a(T0 - T ). Р I dt dx J dx2 F 10 ’

Приведение последнего уравнения к безразмерному виду дает следующее:

ду + ду д2 у

— + — = а—— - ку . дt дх дх

Безразмерные величины в последнем уравнении суть следующее:

X = х, Г = t / (к) ; у = (Т - 70)7(71 - ТО); а = —^-; к = Пак

L \ u J рс • L • u F •peu

Параметр a отражает соотношение интенсивности процессов теплопроводности и переноса тепла движущейся средой, параметр к отражает соотношение интенсивности теплоотдачи в окружающую среду через стенки канала и переносом тепла по каналу.

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

y(t, x = 0) = 0, y(t = 0, x) = 1.

Это соответствует тому, что в начальный момент времени теплоноситель имеет температуру, равную Ti. В момент времени t = 0 в канал начинает поступать жидкость с температурой T0 , что является причиной возникновения переходного процесса. По окончании переходного процесса в системе установится новое равновесие. Целью моделирования является анализ развития переходного процесса.

Представленная выше начально-краевая задача имеет второй порядок по пространственной координате. На входе в канал ( x = 0) краевое условие имеет простой физический смысл: жидкость, поступающая в канал, имеет постоянную температуру. Второе краевое условие в некотором смысле является искусственным. Полагаем, что при x = 1, т.е. на выходе из канала и далее в ближайшей окрестности влиянием теплопроводности можно пренебречь. Тогда краевое условие будет иметь вид

ду ду

-4 + — = -ку.

dt дх

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

Уи+1 Уп Уи+1 Уп+1 Уп+1 2 Уп+1 + Уи+1

yj - yj , yj - yj-1 = ayj+1 - 2 yj + yj -1 ,yn+1

1 : =a ô У i ,

a i 1 Ln+1 I 1 i 1 i 2a , / Ln+1 , a ,.n+1

, A n-\-L , A , 7 П i-1 , ^ П i-1

r„2 + h lyj-1 -\~z+i:+—+к Iyj +—yj+1

''y"

J

j

T

X

Краевые условия:

уп+1 - уп уп+1 - п+1

1 = 1, у'"' = 0. у = М, уу 21 =-ку"+|.

^ Т п ^

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

уМ = х2 уМ-1 +и2 ; 1 = М, у'+1 = 1—п-у"+1 + 1-1-у'' .

— + — + к — + — + к

т к т к

Путем сравнения последних соотношений получим значения коэффи-

циентов Х2 и ^2. Используя основное соотношение прогонки у 1 = а^у 1+1 + в у , получим выражение для ум :

,Л+1 = У2 + Х2вт-1

уМ =~,------------,

1 -Х2 аМ -1

которое используется для расчета остальных значений у1+1. Дальнейшее численное решение задачи проводится по схеме, представленной выше (рис. 6).

Следующей задачей с распределенными параметрами является задача моделирования волновых процессов. В качестве примера рассматривается задача гидроудара. Пусть длинный трубопровод постоянного сечения подсоединен к достаточно большому резервуару с жидкостью. На другом конце трубопровода находится клапан (рис. 8). Под действием перепада давления по трубопроводу течет жидкость с известной скоростью У0. В момент времени t = 0 клапан моментально закрывается. Так как жидкость является упругой средой, то в трубопроводе развивается волновой процесс. Результатом этого является многократное повышение давления в трубопроводе. Предположим, что давление жидкости в резервуаре остается постоянным.

Рис. 8 Схема объекта исследования при гидроударе

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

Краевые условия данной задачи имеют вид

У(ґ, х = Ц) = 0, Р(ї, х = 0) = Ро.

Начальные условия: в момент времени ї = 0 по всему каналу скорость движения жидкости и давление считаем постоянными:

Р(ї = 0, х) = Ро; У(ї = 0, х) = Уо.

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

ЗУ ЗР дР 2 ЗУ п

Р— = —, — + рс — = 0,

дt дх дt дх

где р - плотность жидкости; с - скорость распространения упругих возмущений в жидкости (скорость звука).

Задача имеет пять параметров, поэтому предварительно преобразуем ее к безразмерному виду. В качестве масштаба для продольной координаты выберем длину канала, а в качестве масштаба скорости - скорость звука. В соответствии с этим введем безразмерную координату и скорость: х = х/Ь , У = У / с. В качестве масштаба времени выберем время распространения звука по каналу Ь. Тогда безразмерное время t = t /1 Ь |. С учетом этих соотно-с I с )

шений исходные уравнения примут вид

1 зр зу л зу 1 зр „

+ ^ = 0, ^ + —-^ = 0.

рс2 дt Зх дt рс2 Зх

Для данной задачи величина абсолютного давления в канале Р0 не имеет значения, больший интерес представляют отклонения давления от стационарного значения Р0 :

АР = Р - Р0;

зу + 1 З( р-Рз) = 0- 1 З(Р - Р0) + зу = 0,

дt рс2 Зх ’ рс2 дt Зх ’

ЗАР ЗУ Л ЗУ ЗАР Л

—^ + ^ = 0; —— + —^ = 0,

дt З х дt З х

— 2

где АР = АР / р с . Тогда начальные условия задачи примут вид

аР(1 = 0, х) = 0; У (? = 0, х) = У° = У0.

с

Начальные условия содержат параметр У) - безразмерную начальную скорость движения жидкости. Этот параметр можно исключить, если провести следующие преобразования переменных: У = У / У0; А Р = А Р / У0.

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

ЗУ ЗАР ЗАР ЗУ п т п =-- 1Ч п

—^ + —^ = 0; —-=- + -= = 0. А Р(1 , х = 0) = 0, У (!, х = 1) = 0,

31 З х 31 З х

аР(1 = 0, х) = 0, У (? = 0, х) = 1.

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

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

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

Л/2 -У„1/2 , Ар+1 -ЛР„ = о -Лр+1 + УП+1 -у„ = 0;

т И ’ т И ’

Уу+1/2 = 0,5(У/ + У]+1) ЛРу+1/2 = 0,5(ЛРП +ЛР]+1);

Х]+т = Ху + И/2, х-т = Ху - И/2, Ху+1 = Ху + И.

2. Определение уточненного решения для момента времени ¿п+1 по схеме:

У/+1 -Уу + ЛРу+11/2-ЛР5”+11/2 = 0; ЛР”+1 -ЛР„ + У„:1,2 -Vу%2 = о т и ; т и ,

причем у = 2 (М -1), где М - количество точек по оси х.

В концевых точках канала решение рассчитывается с учетом краевых

условий. Для точки х = 0 (у = 1) расчетная формула построена на основе со-

отношений

„+1 у„+1 - у„ ЛР„+11/2 -лр„+1 лр„+1 =0, ^+—у+^----------у— = 0.

у т И/2

В сечении х = 1 (у = М) расчетная формула строится аналогично:

А1 лр„+1 -лр„ у„+1 - V„-11//22

У„+1 =0, ------}— + ^------у-1/2- = 0.

у т И/2

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

1) расчет V „+1, ЛР „+1 во всех точках с дробными индексами;

2) расчет ЛР„+1, У„+1 для всех внутренних точек с индексами от у = 2 до („ - 1);

3) расчет ЛР„+1, У„+1 в концевых точках: у = 1, у = М.

Индексная форма записи переменных означает следующее:

ЛР/+1 = ЛР(Г„+1, Ху), ЛР„ = ЛР(^„, Ху), ЛРу+1 = ЛР(Г„, Ху+1),

у«+1 = У(^„+ьХу ) , У„ = У^„, Ху ) ; у = Ху ± И , ¿„+1 = ¿„ + т .

Численное решение задачи следует проводить средствами систем программирования, т.к. схема расчета требует выполнения следующего условия: т< И /2 . Результаты моделирования представлены на рис. 9.

0,5 1 1,5 2 2,5 3 3,5 4 4,5 5 5,5 Б 6,5 7 7,5 8 8,5 9 9,5 10

Рис. 9 Результаты моделирования гидроудара: I - A P(t, х = 1), II - V (t, x = 0)

Рассмотренные примеры дают достаточно полное представление об особенностях моделирования процессов с распределенными параметрами.

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

1. Введение в математическое моделирование / под ред. П. В. Трусова. - М. : Логос, 2004. - 440 с.

2. Самарский, А. А. Математическое моделирование: Идеи. Методы. Примеры /

A. А. Самарский, А. П. Михайлов. - М. : Физматлит, 2001. - 320 с.

3. Пасконов, В. М. Численное моделирование процессов тепло- и массообмена /

B. М. Пасконов, В. И. Полежаев, Л. А. Чудов. - М. : Наука, 1984. - 286 с.

4. Фарлоу, С. Уравнения с частными производными для научных работников и инженеров / С. Фарлоу. - М. : Мир, 1985. - 384 с.

5. Королев, А. Л. Аспекты компьютерного математического моделирования / А. Л. Королев // Вузовское преподавание: проблемы и перспективы : материалы по итогам работы 8-й Международной научно-практической конференции (Челябинск, 30-31 октября 2007 г.). - Челябинск : РЕКПОЛ, 2007. - С. 88-96.