Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена.djvu), страница 36

Вычисление РкзультиРувних пОтОИОВ ОО 8 1=1,Н с суимиРовлник поглон!Нных потоков ОО 8 1-[,И 1РП-Л6,7,7 6 К Л«(1-1)/2»1 РУ! Р(К)»3(1)/3(.!) Ркс. 6.3 в)то а 7 Е !»(1-1)/2»» Рз! Р(К) а РНЕ!(1) И(ЕЕ(1)»ЕР$(1)»РЩ»РЕР(Л) веяли ЮВ Рис. 6.3. Продолжение В одномерный массив Е эта информация записывается по столбцам, т. е. в следующем порядке: <рп»р»е <рею <р!» ", (рл — 1,7, (рл„. При этом нетрудно убедиться, что порядковый номер й элемента в одномерном массиве Е, соответствующего коэффициенту срзь рассчитывается по индексам 1,1 так: й = ((( — !)/2 + 1, где 1' ( 1(см. операторы 31, 35, 40).

Если нужно найти коэффициент (рм прн 1 ) (, то сначала выбирается коэффициент (р)1 из массива Г, а затем проводится расчет по формуле (рп = (р(15;/51 (оператор 36). Система уравнений (6.6) решается в подпрограмме с помощью стандартной процедуры СЕ).6.

При этом вектор правых частей В, имеющий смысл собственных потоков, записывается в массив РЕР, в котором после решения системы (на выходе нз С!Е(.б) размещается результат решения — вектор эффективных потоков. Параметры и логическая структура подпрограммы КАП доста. точно подробно пояснены в комментариях к тексту. Для проведения расчетов головная программа, реализующая ввод данных и печать результатов, должна быть составлена пользователем.

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

Система нестационарных уравнений теплового баланса для определения среднеобъемных температур Т; с учетом лучистого теплообмена имеет вид С; — ' =Р! — ~~ о( (Т( — Т)) т Р,'.", 1=-1,, У, (6.9) дт /= 1 где оы — тепловые проводимости, учитывающие конвективный и кондуктивный перенос теплоты; Р; — мощности источников теплоты, действующих в телах; РР(" (Т„..., Тя) — результирующие потоки излучения, зависящие от всех температур н определяемые из решения системы (6.6) и формул (6.8). При использовании для численного решения задачи (6.9) какой- либо явной схемы, например схемы Рунге-Кутта, необходимо для 16! вычисления Р!"' (Т„..., Ти) на каждом временном шаге решать систему алгебраических уравнений (6.6), в которой собственные лучистые тепловые потоки рассчитываются по известным значениям температур с поедыдущего временного слоя.

В случае, когда коэффициенты е; не зависят от температуры, коэффициенты матрицы А постоянны для любого временного шага, а изменяются от шага к шагу только векторы свободных членов В, см. (6.7). Поэтому прн реализации численной схемы с целью экономии машинного времени целесообразно не решать на каждом шаге по времени систему (6.6), а один раз перед началом цикла по времени вычислить обратную матрицу Р = А-', элементы которой обозначим через т(ы.

Тогда выражения для эффективных потоков примут вид и Р'Ф = 0В или Р э= ~Р дт(еточ Т~4 51) /=! а для результирующих потоков Ь=~ =: е, ~ т~~~ тр11 ~ ',~ т(ть еь ач Т4 5ь — о, Т,' 51 (6.10) 1.7=1 ~ь=~ Обратная матрица 0 может быть найдена, например, с использованием стандартной подпрограммы М1)ЧУ (см. 2 !.3).Формирование матрицы А реализовано в приведенной выше подпрограмме (операторы 26 — 43). При использовании для решения системы (6.9) стандартной программы К КО5, реализующей метод Рунге — Кутта четвертого порядка (см. 2 1.6), вычисление правых частей, в том числе расчет Р7' согласно (6.!О), должно быть реализовано в составленной пользователем подпрограмме.

й ь.х Рдсчет хгловых коэефициентов. ВЫЧИСЛЕНИЕ КРАТНЫХ ИНТЕГРАЛОВ Угловой коэффициент излучения 1301 между поверхностями 5, и 5, вычисляется по формуле (6.1 1) з,з где г — расстояние между элементарными площадками Ю, и д5з; 6„0 — угол между лучом, соединяющим площадки, и нормалями к ним (рис. 6.4).

Таким образом, в общем случае для определения угловых коэффициентов необходимо вычислять четырехкратные интегралы (6.11). Для некоторых видов «взаимодействующих» поверхностей интегрирование можно выполнить аналитически. Получающиеся при этом 182 формулы для угловых коэффициентов приведены в литературе [8,30). Однако для поверхности достаточно сложной формы аналитическое интегрирование часто оказывается невозможным, поэтому прибегают к численному интегрированию.

В главе 2 мы познакомились с методами вычисления одномерных интегралов. Г!ри переходе к кратным интегралам возникают новые проблемы, связанные с разбиением области интегрирования и У; зю Рис. 6Л Рис. 6.6 выбором численного метода. Рассмотрим основные методы вычисления кратных интегралов на примере задач расчета угловых коэффициентов. Область интегрирования для (6.11) такова, что ее разбиение на подобласти можно организовать путем независимого дробления двумерных областей, соответствующих поверхностям 5, и 5,.

Поэтому ниже основные методы численного интегрирования будут иллюстрироваться на примере двумерных интегралов. Метод ячеек. Рассмотрим двукратный интеграл ! =) /(х, у)йхс(у (6. 12) в области ьз, представленной на рис. 6.5. При использовании метода ячеек область й разбивается на элементарные подобласти И; и приближенное значение интеграла вычисляется по формуле М 1= ~ ~(хь у~) 5» (6.13) с=! где 5, — плошадь 1-го элемента разбиения Йб х,, у; — координаты точки, принадлежащей 1-му элементу; М вЂ” общее число элементов. Лля использования формулы (6.13) необходимо выбрать вид элементов разбиения и положение точки, в которой вычисляется газ значение интегрируемой функции, внутри элемента.

В качестве элементов обычно используют прямоугольники и треугольники, а точку берут в центре тяжести элемента. При таком подходе формула (6.!3) будет точна для линейной функции ~ (х, у), если область ьз точно разбивается соответствующими элементами. В случае произвольной дважды непрерывно дифференцируемой функции Г (х, у) при одинаковом числе и отрезков разбиения по обеим координатам формула (6.13) имеет второй порядок точности, т.

е. ее погрешность убывает пропорционально 1(л'. бl а Ь а Х Рис. бд На практике область (з часто имеет криволинейную границу, и ее не удается точно разбить на треугольные или четырехугольные элементы. При этом обычно используются достаточно грубые приближения для граничных слагаемых, и тогда в целом погрешность формулы (6.13) имеет порядок О (1ьч). Если подынтегральная функция имеет лишь первые кусочно-непрерывные производные, то формула (6.13) имеет также погрешность порядка О (1!и). Сложная форма области и отсутствие непрерывных вторых производных часто встречаются при расчете угловых коэффициентов.

Для иллюстрации менее очевидного второго факта рассмотрим расчет срлв — со для плоского случая между кривыми АВ и Сх1, изображенными на рис. 6.6, а. Этот угловой коэффициент можно представить как интеграл по прямоугольной области, показанной на рис. 6.6, б: %лв — со= ~)(х, у)дну, где функция ~(х, у) определяется путем проведения соответствующей замены переменных. Остановимся на особенностях этой функции.

Зафиксируем какое-либо значение х, например х = а, и рас- ~вз смотрим поведение функции / (а, у) в зависимости от у. Из рис. 6.6, б видно, что на интервале (с, е! значения /(а, у) будут отличны от нуля, а на интервале !е, ь(! — равны нулю, так как часть кривой ЕР «не видна» из точки А первой кривой. Поскольку при подходе к точке у = е со стороны меньших значений / (а, у) будет убывать до нуля по косннусоидалъному закону, то в этой точке первая производная по у будет иметь разрыв. Рассмотренный пример имеет характерную особенность, обусловлнвающую некоторые сложности, возникающие при составлении программ вычисления угловых коэффициентов.

Она заключается в наличии у поверхностей взаимной затененности, которая может быть вызвана кривизной этих поверхностей или наличием в системе других поверхностей. В последнем случае сама подынтегральная функция в (6.11) может иметь разрывы. Из-за затененности для каждой пары элементов разбиения, принадлежащих разным поверхностям, приходится проверять, «видят» ли они друг друга. В многомерном случае реализация такой проверки прн наличии многих поверхностей может быть достаточно громоздка.

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