Патанкар С. Численные методы решения задач теплообмена и динамики жидкости (1185911), страница 11
Текст из файла (страница 11)
Остальные узловые точки назовем внутренними, вокруг каждой из них размещается контрольный объем. Дискретный аналог, подобный уравнению (3.2), можно записать для каждого такого контрольного объема. Если (3.2) рассматривать как уравнение для определения Тр, то имеем необходимые уравнения для всех неизвестных температур во внутренних узловых точках. Однако два из этих уравнений включают значения температур в граничных узловых точках. Именно через эти температуры на границе осуществляется учет заданных граничных условий в схеме численного решения.
Поскольку нет необходимости обсуждать каждую из граничных точек отдельно, рассмотрим расположенную слева граничную точку В, которая примыкает к первой внутренней точке 1, как показано на рис. 3.3. Обычно в теории теплопроводности на границе могут быть заданы; 1) температура, 2) тепловой поток и 3) тепловой поток, определенный через коэффициент теплоотдачи и температуру окружающей жидкости.
Если задана температура на границе (т. е. значение Тв известно), то не возникает особенных трудностей и не требуются дополнительные уравнения. Когда температура на границе неиз- 44 дг(тв — тг) Чв — (б„) + (Бс+ БРТВ) йх = О. (3.16) Дальнейшее преобразование этого уравнения зависит от того, как задан тепловой поток г)в на границе. Если задано само значение дв, уравнение для Тп записывается следующим образом: а„Тв = агТг + Ь, (3.17) еде Аг аг = (бк); Ь = Всйх + г)в' ав — — аз —:з'рЛх. (3.18) Если тепловой поток определен через коэффициент теплоотдачи Й' и температуру окружающей жидкости Ть такую, что (3.19) то уравнение для Тв записывается следующим образом: авТв = агТ, + Ь, (3.20) где Кг аг = (бк) г Ь = Всйх+ ЬТ,; ав — — аг — Врйх + Ь.
(3.21) ' Напомним, что символ а использовался в гл, 1 для обозначения энтальнпги. Однако его не надо путать с аналогично обозначенным коэффициентом .геплоотдачи. вестна, необходимо дополнительное уравнение для определения Тн. Это уравнение можно получить с помощью интегрирования дифференциального уравнения по половинному контрольному объему, который, как это показано на рис. З.З, примыкает к границе (этот объем лежит только с одной стороны узловой точки В, поэтому к нему относится половина контрольного объема). В увеличенном масштабе этот объем показан на рис. 3.4. Интсгрируя (3.1) по половинному контрольному объему и определяя тепловой поток г(=МТ/г(х, получаем Яв (дк)г г)в г)з+ (5с+ ЯрТв) Лх = О, (3.15) дл где источниковый член линеаризован обычным образом. Тепловой поток д, на поверхности Р 34 Рис. 3.4. Половин- контрольного объема можно записать по ти ный ионтрольныя пу (3.7).
В результате имеем объем вблизи границы Таким же образом можно получить необходимое число уравнений для неизвестных температур. Далее рассмотрим метод решения этих уравнений. Решение линейных алгебраических уравнений. Решение дискретного аналога для одномерного случая можно получить с помощью стандартного метода исключения Гаусса.
Для уравнений такого простого вида процесс исключения превращается в очень удобный алгоритм. Иногда оп называется алгоритмом Томаса илн ТВМА (Тг)-с)!азиопа(-Ма1пх А!дог!йп> — трсхдиагональный матрицы алгоритмом) '. Название ТРМА является результатом того, что, когда матрица коэффициентов этих уравнений записана, все ненулевые коэффициенты группируются вдоль трех диагноналей матрицы.
Для удобства записи алгоритма введем некоторые обозначения. Присвоим узловым точкам, изображенным на рис. 3.3, номера 1, 2, 3, ..., Х Номера ! и !»' относятся к точкам на границе. Дискретный аналог можно записать в следующем виде: а>Т, = Ь,,Тг»И+ с,Тт, + с(п (3.22> где >=1, 2, 3, ..., >тг Таким образом, температура Те связана с соседними значениями Тесы и Те ь Запись уравнений для узловых. точек на границе лает п>=О и Ьм=О, следовательно, температуры Т, и Тмч» не будут иметь смысла (в том случае, когда температуры на границе заданы, уравнения для граничных точек записываются в обычной форме, например если Т, задано, имеем а,=1; Ь>=0, с,=О и с1, равно заданному- значению Т>). Записанные условия означают, что Т, известна в зависимости от Та.
Уравнение для 1=2 представляет собой соотношение между Ть Т, и Т,, Но поскольку Т> может быть выражена через Т>, это соотношение приводится к соотно>пени~о между Т, н Т,. Другими словами, Т, можно выразить через Тз. Процесс подстановки можно продолжать до тех пор, пока значение Т» не будет выражено через Т, еь Но поскольку Тат > ие существует, мы в действительности на данном этапе получим численное значение Тм. Это позволит нам начать процесс обратной подстановки, в котором Т.,=, получится из Та>, Та=а — из Т» ь ..., Те — из Т, и Т> — из Т,.
Это и составляет существо алгоритма трехдиагоиальной матрицы. Предположим, что при прямой подстановке имеем зависимость Те = Р,.те ы + (;>е (3.24) после того, как получено Tг — > =- Рг->Т;+ !',>е — >. (3.25) ' В отечественной технической литературе этот метод носит название метода прогонки. — Прим. науч. ред. Подставляя (3.25) в (3.22), получаем следующее соотношение: а,.Т,. = Ь,.Тс ы + с,. (Рс, Т, + Я; 1) + с(п (3.26) которое можно привести к виду (3.24).
Иначе говоря, коэффициенты Р; и Яс запишем в виде ь; Р,= гч — с;Р; 1 (3.27) тз е;+ сдс— с а; — с;Р; .Эти рекуррентные соотношения определяют Р; и Яс через Р;, и Я; ь Заметим, что в начале рекуррентного процесса уравнение (3.22) для 1=1 по форме почти совпадает с (3.24). Таким образом, Р, и Я, определяются в следующем виде: Р, = Ь,/а,; 1г, = с(,/а,. (3.28) (3.29) С этого момента осуществляется обратная подстановка с помощью уравнения (3.24). Краткое описание алгоритма 1. Рассчитываем Р, и Я~ нз уравнения (3. 28). 2.
Используя рекуррентные соотношения (3.27), получаем Р; и Яс для 1=2, 3, ..., Лс. 3. Полагаем Тн — — гн. 4. Используя уравнение (3.24) для с=Лс — 1, Лс — 2, ..., 3„2, 1, по.лучаем Тн ь Тн-ж, Тз, Тм Ть Алгоритм, использующий свойство трехдиагональности матрицы, является мощным и удобным методом решения алгебраических уравнений„которые можно представить в виде (3.22).
В отличие от общих матричных методов Т11МА требует машинной памяти и машинного времени, пропорциональных У, а не Л1' или Л1з. з.з. нестАцнОнАРнАВ ОднОмеРнАя теппОпеОВОднОсть Обобщенный дискретный аналог. Что касается общего дифференциального уравнения для Ф, то теперь мы умеем, по крайней мере, в одномерном случае аппроксимировать диффузионный и источниковый члены. Обратимся к нестационарному члену и временно опустим источниковый член, так как в данном случае нет [Интересно заметить, что эти выражения следуют из уравнения (3.27) после подстановки С1 =0.] На другом конце последовательности Ро Яс имеем Ьн=О.
Это .дает Рн=О, и из (3.24) получаем Тн = Ян. необходимости говорить о нем, Таким образом, ищем решение нестационарного одномерного уравнения теплопроводности (3.30) В дальнейшем для удобства будем полагать рс постоянным (в гл. 1 показано, как уравнение теплопроводности можно модифицировать с учетом переменной теплосмкости). Поскольку время является однонаправленной координатой, решение получаем, передвигаясь во времени от заданного начального распределения температуры. Таким образом, на типичном временнбм шаге по заданным значениям Т в узловых точках для времени с надо определить значения Т для времени г+сзс.
Старые (заданные) значения Т в узловых точках обозначим Тр', Тв~, ТсР„ а новые (неизвестные) значениЯ длЯ вРемени с+ЛС вЂ” Тг', Тв', Ти-'. Дискретный аналог получим путем интегрирования уравнения (3.30) по контрольному объему, показанному на рис. 3.2, и по времсннбму интервалу от с до С+Ы. Таким образом, е с+ы с+ы 1 рс ~ ~ — с((с1х = ~ ~ — (й — ) с(х Ю, (3.31) где пределы интегрирования выбраны в соответствии с физическим смыслом членов. Для представления члена дТ)д( предположим, что значение Т в узловой точке распространено на весь контрольный объем, тогда е с+ы рс ~ ~ — Йс(х = рсбх(Тр — Тр) . бс (3.32~ Следуя способу аппроксимации члена й дТ/дх в стационарном слу- чае, получаем с+ы рсйх (Т~ — Т~~) = ( ~ ( " 1 Ш.
(3.33) (бх) (бх),„ На данном этапе необходимо ввести предположение относительно изменения во времени от с до с+И температур Тн, Тх и Тч. Возможны различные предположения, и одно из пих имеет следующий вид: с+ы ~тТ В +)- — ~))Т',1, (3,34)с где )' — весовой коэффициент, изменяющийся от 0 до 1, Использусь аналогичные соотношения для интегралов от Тх до Тж из урав- нения (З.ЗЗ), находим 48 ( йе (т,'.— тр] й (т„'— тД р» — (Тр — Тр/ = / ~ д) ~ (бх), (бх) . ~+ + (1 — /) ' йе(т» — т',) й (тр — т'] (бх)е (бх)м (3.35) Преобразуя это выражение, опустим индекс 1'и запомним, что Тр, Тн, Т,г с этого момента будут означать новые значения Т для времени /+Лб В результате имеем а„ТР = ав ')'/Т» + (1 — /) То] + а, (/Тн + (1 — /) То ] + + (ао — (1 — /)а.— (1 — /)а,] Тр, (З.3бт где а . =)((,/(бх),; а, = й„,/(бх); ар =Рсбх/М; ар =/а +/а, + ао.
(3.37) а Тр = а Т', + а„,Т'„, + (а' — а — а, )Т'. (3.38) Это означает, что Тр не зависит от других неизвестных, таких, как Тк или Т„,, а является явно определенной по известным тем- Явная, Кранка — Николсона н полностью неявная схемы. Для определенных конкретных значений весового коэффициента дискретный аналог приводится к хорошо известным схемам для параболических диф- ", г 2 ференциальных уравнений.
В частности, для тр 1 1=0 получаем явную схему, для /=0,5— схему Кранка — Николсона и для /=1 —. полностью неявную схему. Кратко рассмотрим эти схемы и покажем, что неявная схема наиболее предпочтительна. Различные значения / можно интерпре- и тировать как характеристику изменения Тр от Г, показанного на рис.
3.5. Явная схема Рнс. Зд. Изменение темпо существу предполагает, что старое зна- "ературы но еремеем чение Тр сУЩествУет в пРеделах всего схемы Кранкн — Инке. со- о временнбго шага, за исключением точки нн (й) й полностью не- /+Ы. Неявная схема предполагает, что явной схемы (е) в момент г Тв резко уменьшается от Тр' до Тр', а затем остается равной Тг' на всем временнбм, шаге и температура в пределах временнбго шага характеризуется новым значением Тр', Схема Кранка — Николсона предполагает линейное изменение Тр, С первого взгляда линейное изменение.
должно быть более разумным, чем две другие альтернативы. Почему же мы предпочитаем неявную схемуй Ответ будет очень коротким. Для явной схемы (/=О) уравнение (3.36) принимает следующий вид: -леРатУРам Тв", Тка, Т„'. ПозтомУ схема и называетсЯ Явной, Любая схема с /~=0 должна быть неявной, так как Тг зависит от неизвестных Тв и Т„,, в этом случае необходимо решать одновременно несколько уравнений. Удобство явной схемы в этом отношении компенсируется, однако, рядом ограничений, Лналпзируя (3.38) и вспоминая основное правило о положительных коэффициентах (правило 2), замечаем, что коэффициент при Тг' может принимать отрицательные значения (значение Тг' рассматривается как соседнее с Тг по временнбй координате).