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

Известно, что решения задач теплопроводности, полученные с помощью классических аналитических методов, представляются в форме бесконечных рядов, плохо сходящихся в окрестностях граничных точек и при малых значениях временной координаты. Исследования показывают, что сходимость точного аналитического решения нестационарной задачи теплопроводности для бесконечной пластины при граничных условиях первого рода в диапазоне чисел Фурье 10 12 < Fo < 10 7 наблюдается лишь при использовании от 1000 (Fo = 10 7) до пятисот тысяч (Fo =10 ,2) членов ряда.

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

В аналитической теории теплопроводности известны методы, в которых используется понятие глубины термического слоя (интегральные методы теплового баланса) [8, 9, 13, 19, 38, 67, 86|. При их использовании процесс нагрева (охлаждения) тел формально разделяется на две стадии. Первая из них характеризуется постепенным продвижением фронта температурного возмущения от поверхности к центру, а вторая — изменением температуры по всему объему тела вплоть до наступления стационарного режима. Перемещение фронта температурного возмущения учитывается введением новой функции q(Fo)> называемой глубиной проникания (глубиной термического слоя). Такая модель процесса теплопроводности используется в ряде методов: интегральном методе теплового баланса [8,19]; методе осреднения функциональных поправок [67]; методах Швеца [86], Био [9], Вейника [13] и др.

Несомненным их преимуществом является возможность получения простых но форме аналитических решений удовлетворительной точности как для регулярного, так и нерегулярного процессов теплопроводности. Однако их серьезным недостатком является низкая точность. Причина в том, что получаемое решение, точно удовлетворяя начальному и граничным условиям, основному дифференциальному уравнению удовлетворяет лишь в среднем. Это связано с тем, что в основу метода положено построение так называемого интеграла теплового баланса, что равнозначно осреднению исходного дифференциального уравнения в пределах глубины термического слоя. Следовательно, очевидным путем повышения точности интегральных методов является улучшение выполнения исходного дифференциального уравнения. С этой целью в настоящей работе избрано направление аппроксимации температурной функции полиномами более высоких степеней. Для определения неизвестных коэффициентов таких полиномов основных граничных условий оказывается недостаточно. В связи с этим возникает необходимость привлечения дополнительных граничных условий, определяемых из исходного дифференциального уравнения с использованием основных граничных условий и условий, задаваемых на фронте температурного возмущения. Таким путем можно получать аналитические решения ряда краевых задач практически с заданной степенью точности во всем диапазоне времени нестационарного процесса (0 < Fo < °о) без каких-либо ограничений на величину числа Фурье в области малых его значений.

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

где 0 = (Г - 7о)/(Гст - Го) — относительная избыточная температура; Fo = = ax/R2 число Фурье; ?, = x/R — безразмерная координата; х — координата; R — половина толщины пластины; а — коэффициент температуропроводности; т — время; Г0 — начальная температура; Гст — температура стенки при ? = 0.

Процесс нагрева разделим на две стадии по времени: 0 < Fo < Fo и Fo < Fo < оо. Для этого введем движущуюся во времени границу (фронт температурного возмущения), разделяющую исходную область 0 < ^ < 1 на две подобласти: 0 < ? < ^(Го) и q(Fo) < Ъ, < 1, где q(Fo) — функция, определяющая продвижение границы раздела во времени (рис. 7.9). При этом в области, расположенной за фронтом температурного возмущения, сохраняется начальная температура. Первая стадия процесса заканчивается при достижении движущейся границей центра пластины (? = 1), т.е. когда Fo = Fo. Во второй стадии изменение температуры происходит по всему объему тела 0 < ? < 1. Здесь вводится в рассмотрение дополнительная искомая функция ^2(Го) = 0(1, Го), характеризующая изменение температуры во времени в центре пластины.

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

Расчетная схема теплообмена

Рис. 7.9. Расчетная схема теплообмена

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

Математическая постановка краевой задачи для первой стадии процесса имеет вид

где соотиошеиия (7.73), (7.74) представляют условия сопряжения прогретой и непрогретой зон. Соотношение (7.73) устанавливает равенство температуры тела в точке ?, = q(Fo) его начальной температуре. Согласно условию (7.74), тепловой поток не распространяется за пределы фронта температурного возмущения (условие адиабатной стенки). Математическое доказательство условий (7.73), (7.74) дано в [31 ].

В связи с принятием допущения о равенстве температуры тела на фронте температурного возмущения ?, = q{Fo) его начальной температуре, обратим внимание на тот очевидный факт, что на первой стадии процесса задача (7.71)—(7.74) за пределами фронта температурного возмущения, т.е. на отрезке q(Fo) < ?, < 1, вообще не определена. В связи с этим здесь нет необходимости выполнения начального условия вида (7.68) по всей толщине пластины (поэтому такое условие отсутствует в математической постановке задачи (7.71)—(7.74)). В данном случае вполне достаточным является выполнение граничного условия (7.73), согласно которому для всех ^ = q(Fo) температура тела равна начальной температуре. Кроме того, в данной задаче отсутствует также граничное условие вида (7.70), ввиду того что оно не влияет на процесс теплообмена в первой его стадии.

При использовании классического аналитического решения, полученного, например, с использованием метода разделения переменных, наибольшие проблемы возникают в случае нахождения температуры для малых значений числа Фурье (Fo5? 0). Это связано с тем, что искомое решение в данном случае должно описывать две линии — прямую DC и линию АЕ, имеющую незначительную кривизну (см. рис. 7.9), связанные между собой криволинейным участком ED и находящиеся почти под прямым углом друг к другу. Причем на одной из этих прямых (линия DC) температура неизменна, а на другой при незначительно изменяющейся координате ?, температура изменяется практически от нуля до единицы. Удовлетворить всем этим столь неоднородным и изменяющимся во времени условиям в одном аналитическом выражении (при Fo —? 0) можно лишь при использовании в нем бесконечно большого числа членов ряда, о чем уже упоминалось выше.

Преимущество данного метода состоит в том, что получаемое решение в данном случае не связано с необходимостью аппроксимации температуры на участках DC, ВС и др., т.е. за фронтом температурного возмущения q(Fo) < ^ < 1, так как в этой области на первой стадии процесса задача (7.71)—(7.74) не определена. Получаемое аналитическое решение описывает изменение температуры во времени, характеризуемое лишь кривыми вида AD, АВ и др., для практически точного описания которых достаточно всего нескольких членов ряда решения вне зависимости от величины врс- меннбй координаты.

Отметим, что задача (7.71)—(7.74) нс относится к классу задач, в которых учитывается конечная скорость продвижения тепловой волны. Получение их решений сводится к интегрированию гиперболического (волнового) уравнения теплопроводности [49]. Введенный в задаче (7.71)—(7.74) фронт температурного возмущения по физическому смыслу является аналогом движущейся изотермы (но не тепловой волны). Ввиду того, что на фронте температурного возмущения в процессе его движения по координате Е, поддерживается начальная температура С-)(<7[, Во), он, следовательно, является аналогом нулевой изотермы (ниже будет показано, что скорость продвижения нулевой изотермы при п —*? оо приближается к бесконечному значению, т.е. Fo —*? 0).

Решение задачи (7.71)—(7.74) разыскивается в виде полинома

где ak{(h ) — неизвестные коэффициенты. Для их определения используются граничные условия (7.72)—(7.74). Подставляя (7.75) (ограничиваясь тремя членами ряда) в (7.72)—(7.74), для определения яД#]) (k = 0,1,2) будем иметь систему трех алгебраических линейных уравнений. После определения ait(q) соотношение (7.75) принимает вид

Для нахождения неизвестной функции q{Fo) в первом приближении составим невязку уравнения (7.71) и проинтегрируем ее в пределах глубины термического слоя (что равнозначно построению интеграла теплового баланса — осреднению уравнения (7.71)):

Подставляя (7.76) в (7.77), после определения интегралов получаем уравнение

Интегрируя уравнение (7.78) при начальном условии q(Q) = 0, получаем уравнение

Положив Fo ~ 0,083333.

Соотношения (7.76), (7.79) определяют решение задачи (7.71)—(7.74) в первом приближении.

Результаты расчетов по формуле (7.76) в сравнении с точным решением [49] приведены на рис. 7.10. Их анализ позволяет заключить, что расхождение с точным решением составляет 3—4%. При этом основная погрешность возникает из-за неточного выполнения дифференциального уравнения (7.71). Отметим, что граничные условия (7.72)—(7.74) и интеграл теплового баланса (7.77) выполняются точно.

Изменение относительной избыточной температуры в пластине

Рис. 7.10. Изменение относительной избыточной температуры в пластине:

------по формуле (7.76); о — по формуле (7.90);--точное решение

Очевидным путем повышения точности решения является увеличение степени аппроксимирующего полинома (7.75). Для определения появляющихся при этом дополнительных неизвестных коэффициентов необходимо привлекать дополнительные граничные условия. Для их получения будем последовательно дифференцировать граничные условия (7.72)—(7.74) по переменной Fo, а уравнение (7.71) — по переменной Сравнивая получающиеся при этом соотношения, можно найти любое необходимое количество дополнительных граничных условий. Для получения первого из них продифференцируем условие (7.72) по переменной То:

Сравнивая (7.80) с уравнением (7.71), находим первое дополнительное граничное условие

Продифференцируем соотношение (7.73) но переменной Fo. Так как из (7.73) требуется находить значение ©(^, Fo) в точке Ъ, = q(Fo)> то ?, является функцией Fo и, следовательно, ©(?,, Fo) будет сложной функцией. Тогда

360

Соотношение (7.82) с учетом ^7.741 ппиволится к ни/iv

Сравнивая (7.83) с уравнением (7.71), применительно к точке ?, = q(Fo) получаем второе дополнительное граничное условие

Продифференцируем соотношение (7.74) по Fo с учетом того, что переменная ?, является функцией Fo:

Соотношение (7.85) с учетом (7.84) приводится к виду

Продифференцируем уравнение (7.71) по переменной ?, и применим полученное соотношение к точке ?, = q(Fo):

Так как 0(^, Fo) и получаемые частные производные являются непрерывными функциями, то правомерна замена порядка дифференцирования, т.е.

Сравнивая (7.86) и (7.87), находим третье дополнительное граничное условие

Во втором приближении, используя дополнительные граничные условия (7.81), (7.84), (7.89) совместно с заданными условиями (7.72)—(7.74), можно найти уже шесть коэффициентов полинома (7.75) и задать температурную функцию в виде полинома пятой степени. Подставляя (7.75), ограничиваясь шестью членами ряда, во все перечисленные условия относительно неизвестных коэффициентов получаем систему шести алгебраических линейных уравнений. Подставляя найденные из решения этой системы коэффициенты (k = 0, 5) в соотношение (7.75), получаем соотношение

361

Подставляя (7.90) в (7.77), относительно q(Fo) приходим к обыкновенному дифференциальному уравнению

Интегрируя (7.91) при начальном условии

Положив в (7.92) q(Fo) = 1, находим Fo = = 0,005.

Соотношения (7.90), (7.92) определяют решение задачи (7.71)—(7.74) с дополнительными граничными условиями (7.81), (7.84), (7.89) во втором приближении. Результаты расчетов по формуле (7.90) в сравнении с точным решением [49] приведены на графиках рис. 7.10. Их анализ позволяет заключить, что в диапазоне чисел Фурье 1 • 10 5 < Fo < Fo = 0,05 отклонение полученного по формуле (7.90) решения от точного составляет 1 —2%.

Применяя рассмотренный выше способ, можно получить какое угодно число дополнительных граничных условий. Например, для получения четвертого дополнительного граничного условия продифференцируем граничное условие (7.81) по переменной Fo:

Продифференцируем уравнение (7.71) соответственно два и три раза по переменной Ь;.

Сравнивая (7.93) и (7.94), получаем четвертое дополнительное граничное условие

Продифференцируем дополнительное граничное условие (7.84) по переменной Fo, считая, что = q(Fo) является функцией Fo:

Последнее соотношение с учетом (7.89) примет вид

Сопоставляя (7.94) и (7.97), применительно к точке ^ = q{Fo) находим пятое дополнительное граничное условие

362

Дифференцируя условие (7.89) но переменной Fo, получаем соотношение Соотношение (7.99) с учетом (7.98) примет вид

Сопоставляя (7.95) и (7.100), находим шестое дополнительное граничное условие

Дополнительные граничные условия (7.81), (7.84), (7.89), (7.96), (7.98), (7.101) совместно с заданными условиями (7.72)—(7.74) позволяют найти девять неизвестных коэффициентов а^ (к = 0, 8) из (7.75). Подставляя найденные из решения соответствующей системы алгебраических линейных уравнений значения коэффициентов в (7.75), получаем соотношение для определения температуры в третьем приближении:

Подставляя (7.102) в (7.77), получаем уравнение

Интегрируя последнее уравнение при начальном условии #i(0) = 0, получаем соотношение

Время окончания первой стадии процесса Fo = 0,03472.

Следующие три дополнительных граничных условия имеют вид

Уравнение для q(Fo) в данном случае будет иметь вид Его решение примет вид

Вычисляя Fo при q(Fo) = 1, находим Fo = 0,0265152.

Соотношение (7.75) в четвертом приближении принимает вид

)

где Л0 = 1; А х = -55/16; А2 = А4 = А6 = 0; Л3 = 165/16; Л5 = -231/8; А1 = 825/8; Л8 = -165; Л9 = 1925/16; Л10 = -44; Аи = 105/16.

Для получения решения в пятом приближении необходимо добавить еще три дополнительных граничных условия:

Уравнение относительно q(Fo) в пятом приближении будет иметь вид Его решение примет вид

Время окончания первой стадии процесса Fo = Fot = 0,0214286.

Решение задачи (7.71)—(7.74) в пятом приближении записывается в виде

где А0 = 1;А2 = А4 = А6 = А8 = 0; Л, = -245/64; Л3 = 455/32; Л5 = -3003/64; Л7 = 2145/16; Л9 = -35035/64; Л10 = 1001; Ли = -28665/32; Л12 = 455; Л,3 = = -8085/64; Л14 = 15.

Аналогично можно записать дополнительные граничные условия и для других приближений. В частности, были также получены решения в десятом и четырнадцатом приближениях. Функция q(Fo) для этих приближе- ний соответственно определяется по формулам q(Fo) = 2V82G2/'«/19; q(Fo) = 2л/287Ео/3. Время окончания первой стадии процесса для десятого п четырнадцатого приближений соответственно будет равно Fo = 0,010992 и Fo = 0,00784.

Формула для безразмерной температуры в десятом приближении имеет вид

где Л0 = 1; Л2 = Л4 = А6 = Л8 = Л10 = Л12 = Л14 = Alfi = Л|8 = 0;

Л4 = -352 495/65 536; Л3 = 1 306 305/32 768;

Л5 = -1 6981 965/65 536; Л7 = 21 460 725/16 384;

Л9 = -350 525 175/65536; Аи = 605 452 575/32 768;

Л13 = -3 732 515 325/65 536; Л15 = 1 386 362 835/8192;

Л17 = -37 105 593 525/65 536; Л19 = 121 732 385 775/32 768;

Л20 = -10 015 005; Л21 = 991 249 427 025/65 536;

Л22 = -15 607 800; Л23 = 191 981 114 325/16 384;

Л24 = -6 531 525; Л25 = 176 622 625 179/65 536;

Л26 = -803 880; Л27 = 5391 411 025/32 768;

Л28 = -20 735; Л29 = 79 676 025/65 536.

Отметим, что соотношения для полученных выше решений представляют степенные алгебраические полиномы относительно переменных с, и Fo, которые не содержат как специальных функций (Бесселя, Лежандра, гамма-функций и др.), так и тригонометрических.

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

Результаты расчетов для третьего, седьмого, и четырнадцатого приближений в сравнении с точным решением [49] даны на рис. 7.11. Их анализ приводит к заключению о том, что с увеличением числа приближений решение всякий раз уточняется. Так, уже в седьмом приближении значения температур в диапазоне 5 10 9 < Fo < Fo отличаются от их точных значений не более чем на 0,002%, а в четырнадцатом приближении — на 0,0004%. Следует отметить трудности получения точного решения по формулам из [49] для столь малых чисел Фурье ввиду необходимости использования большого числа членов ряда точного решения. В частности, расчеты показали, что при Fo = 10 7 для сходимости точного решения необходимо использовать около 1000 членов ряда (см. формулу (16), с. 87 из [49]). Для чисел Fo = 10 8; 10'9; 10 10; 10 п; 10 12 сходимость точного решения наблюдается соответственно при величинах чисел ряда 5000; 10 000; 50 000; 200 000; 500 000.

Изменение относительной избыточной температуры в пластине

Рис. 7.11. Изменение относительной избыточной температуры в пластине:

  • --------— третье приближение;-------седьмое;
  • -------------четырнадцатое;--точное решение

Анализ полученных результатов позволяет сделать заключение о низкой эффективности применяемой в классических методах линейной суперпозиции частных решений с целью выполнения начального условия. Именно на этом этапе получения классического решения происходит его максимальное усложнение, что связано с необходимостью подчинения решения столь неоднородным условиям, о которых уже упоминалось выше и которые оказываются практически невыполнимыми при Fo —> 0 (Fo ^ 0) ввиду необходимости использования в классическом решении бесконечно большого числа членов ряда. Представление исходной краевой задачи в виде двух взаимосвязанных процессов, рассматриваемых раздельно и связанных лишь условием сопряжения при Fo = Fo{, позволяет избежать указанных трудностей при возможности получения решения практически с любой заданной степенью точности.

Повышение точности решения с увеличением числа приближений происходит за счет увеличения точности выполнения уравнения (7.71), что подтверждается анализом изменения его невязки в (рис. 7.12). Отметим, что невязка уравнения (7.71) для любой точки координаты ^ определяется лишь начиная с момента, когда фронт температурного возмущения достигает этой точки, и заканчивая моментом времени окончания первой стадии

Изменение невязки е уравнения (7.71) от безразмерного времени Fo

Рис. 7.12. Изменение невязки е уравнения (7.71) от безразмерного времени Fo

в точке % = 0,5:

о — седьмое приближение (Fo = 0,0154762); ? — четырнадцатое приближение (Fot = 0,00784)

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

С целью дополнительного, более глубокого анализа получаемых с помощью интегрального метода решений проведем исследование закономерности изменения фронта температурного возмущения q(Fo) во времени в сравнении с точным решением [49]. Графики перемещения q(Fo) по координате ?, во времени даны на рис. 7.13. Их анализ позволяет заключить, что с увеличением числа приближений время Fo достижения фронтом температурного возмущения координаты ^ = 1 уменьшается. Так, например, в первом приближении оно составляет Fo = 0,0833. Для данного Fo, согласно точному решению [49], температура в центре пластины (? = 1) будет равна 0Ц(1, Fo = 0,0833) = 0,97139. Согласно решению интегральным методом, температура в данной точке равна 0(1, Fo) = 1 (так как в интегральном методе температура на фронте температурного возмущения, в том числе и в точке ?, = 1 при Fo = Fo, принимается равной начальной температуре 0(^, 0) = 1). Принятие начального условия равным единице, а граничного условия первого рода равным нулю (в отличие от задачи (7.67)—

(7.70)) связано с удобством сравнения с точным аналитическим решением, полученным именно в такой постановке [49]. Очевидно, что разность v = = 0(1, Fo) - ©ц(1, Fo{) = 1 - 0,971 = 0,0248 (2,48%) является отклонением

Кривые перемещения фронта температурного возмущения q{Fo) по координате ? во времени Fo

Рис. 7.13. Кривые перемещения фронта температурного возмущения q{Fo) по координате ? во времени Fo:

1, 2, 3,4, 5, 7, 10, 14 — номера приближения решения, полученного интегральным методом, от точного аналитического решения.

Время окончания первой стадии процесса во втором приближении составляет Fo = 0,05. Для этого значения температура, полученная из точного решения, будет равна 0Ц( 1, Fo) = 0,9969. Следовательно, погрешность второго приближения по отношению к точному решению составит v = 1- 0,9969 = = 0,0031 (0,31%). В третьем приближении расхождение решений будет равно v = 1 - 0,9997 = 0,000295 (0,0295%), в четвертом приближении — v = = 0,28177 • 10 4 (0,0028177%), в четырнадцатом — v = 0,278538 • 10 14 (0,278538-10 12%).

Значения температур 0„(1, /’о,), полученных из точного аналитического решения, для различных значений чисел Fo, найденных интегральным методом, а также величина v, приведены в табл. 7.4.

Таблица 7.4

Номер

приближения

1

2

3

4

5

Fo,

0,08333

0.05

0,0347222

0,0265152

0,0214286

©ц(1, Fo,)

0,97139

0,996869

0,9997044

0,99997182

0,9999972756

V

0,028608

0,003131

0,2956 ? 10 3

О

to

  • 00
  • —1

о

0,2724-10 3

Номер

приближения

7

10

14

Fo,

0,0154762

0,0109195

0,00783972125

0„(1, Fo,)

0,999999973681

0,99999999997367

0,9999999999999972146

V

0,263185416 • 10 7

0,26331566 -10 10

0,278538 • 10 14

Анализ данных табл. 7.4 позволяет заключить, что при неограниченном приближении температуры в центре пластины ©u(l, Fo), определяемой из точного аналитического решения [49], к начальной температуре 0(6,, 0) = 1, т.е. при 0(%, 0) = 0Ц(1, Fo), величина Ft^ приближается к нулевому значению (Fot —*? 0). Полученный результат полностью согласуется с гипотезой о бесконечной скорости распространения теплового возмущения, лежащей в основе вывода параболического уравнения теплопроводности вида (7.67). Согласно этой гипотезе, с момента начала действия граничного условия 0(0, Fo) = 0 при = 0 температура на всем отрезке координаты 0 < Ъ, < 1, в том числе и в центре пластины (6, = 1), уже не равна начальной температуре 0(^, 0) = 1, а отличается от нее на некоторую (а в центре пластины на бесконечно малую) величину А0 (рис. 7.14). Значения этой величины, полученной из расчетов по точному аналитическому решению, в зависимости от координаты (0 < ^ < 1) для Fo = 0,005 даны в табл. 7.5.

Таблица 7.5

6

0

0,1

0,2

0,3

0,4

0,5

Д0 = 1 - 0(6, Fo)

1

0,3173

0.0455

0,0027

0,0000633

0,5733 • НГ6

6

0,6

0,7

0,8

0,9

1

Д0 = 1 - 0(6, Fo)

0,197 -10 8

0,266-10 11

0,124-10 14

0,15 • 10 18

0

Отклонение температуры А© = 1 — (?, Fo) на фронте температурного возмущения и за его пределами от начальной температуры

Рис. 7.14. Отклонение температуры А© = 1 — (?, Fo) на фронте температурного возмущения и за его пределами от начальной температуры

(без соблюдения масштаба).

Примечание: числовые значения Д0 для Fo = 0,005 даны в табл. 7.2.

Анализ результатов табл. 7.5 позволяет заключить, что уже при ^ = 0,3 безразмерная температура отличается от начального условия на величину 1 - 0,9973 = 0,0027. В связи с этим на всех графиках, приведенных в известной литературе, считается, что 0(0,3; 0,005) ~ 1, в том числе и для всех > 0,3[49]. Отметим, что 1 - 0(1; 0,005) —? 0 при = 1.

На рис. 7.14 представлены результаты расчетов, выполненных с использованием точного аналитического решения с целью оценки отклонения температуры по толщине пластины (см. табл. 7.5) от начальной температуры А0 = 0(^, 0) - ©(%; 0,005) для Fo = 0,005, где 0(^, 0) = 1.

Из анализа полученных результатов следует, что с увеличением точности определения температуры в интегральном методе время Fo достижения фронтом температурного возмущения координаты ?, = 1 уменьшается и в пределе при п —? °о Fo —* 0. Следовательно, к нулю будет приближаться и расхождение приближенного и точного решений. Анализ изменения температуры и невязки с уравнения (7.71) для различных приближений интегрального метода подтверждает данное заключение (см. рис. 7.11, 7.12).

Метод дополнительных граничных условий можно применить и для второй стадии процесса нагрева (охлаждения). Вторая стадия теплового процесса, соответствующая времени Fo > Fo, характеризуется изменением температуры уже по всему сечению пластины вплоть до наступления стационарного состояния. Для этой стадии понятие термического слоя теряет смысл, и в качестве дополнительной искомой функции принимается функция 0(1; Fo) = q2(Fo), характеризующая изменение температуры от времени в центре пластины (см. рис. 7.9).

Математическая постановка задачи для второй стадии процесса имеет вид

Задача (7.110)—(7.113) не содержит начального условия, что связано со следующими обстоятельствами. При Fo = Fo q(Fo) = 1 и #2(^°i) = 0. Граничные условия (7.111)—(7.113) в этом случае становятся идентичными граничным условиям (7.72)—(7.74), и, следовательно, математические постановки задач (7.71)—(7.74) и (7.110)—(7.113) полностью совпадают. Таким образом происходит плавный переход от первой стадии процесса ко второй.

Так как при Fo = Fo qx(Fo) = 1, то соотношение (7.76) будет иметь вид

Соотношение (7.114) представляет начальное условие задачи (7.110)—

(7.113). Однако в его специальном удовлетворении нет никакой необходимости, так как оно будет выполнено в процессе получения решения этой задачи.

Как и в первой стадии, решение задачи (7.110)—(7.113) разыскивается в виде полинома п-й степени

Неизвестные коэффициенты (k = 0, 1, 2) находятся из граничных условий (7.111)—(7.113). После их определения и подстановки в (7.115) получаем соотношение

Для нахождения решения в первом приближении составим невязку дифференциального уравнения (7.110) и проинтегрируем ее в пределах от ^ = 0 до ^ = 1, т.е.

Подставляя (7.116) в (7.117) и определяя интегралы, относительно неизвестной функции <72(^0)11Рих°Дим к обыкновенному дифференциальному соотношению

Интегрируя уравнение (7.118) при начальном условии q2(Fo) = 0, получаем соотношение

Подставляя (7.119) в (7.116), получаем соотношение

где Fo = 0,0833 (найдено в первом приближении первой стадии процесса).

Результаты расчетов по формуле (7.120) в сравнении с точным решением [49] даны на рис. 7.15. Их анализ позволяет заключить, что максимальное отличие температур, полученных по формуле (7.120), от их точных значений составляет 8%. Анализ решения (7.120) показывает, что оно точно

О 0,2 0,4 0,6 0,8 4 1,0

Рис. 7.15. Изменение относительной избыточной температуры во второй стадии процесса:

--первое приближение;- — второе приближение; Д — точное решение

удовлетворяет начальному условию (7.114) и граничным условиям (7.111)—(7.113). Следовательно, основная погрешность решения происходит от неточного выполнения дифференциального уравнения (7.110). В самом деле, как это следует из соотношения (7.117), уравнение (7.110) удовлетворяется лишь в среднем по толщине пластины.

Увеличение точности решения связано с увеличением числа членов ряда (7.115). Появляющиеся при этом дополнительные неизвестные коэффициенты bk могут быть найдены из дополнительных граничных условий. Для их определения продифференцируем граничные условия (7.111)—(7.113) по переменной Fo:

Сравнивая соотношения (7.121) и (7.122) с уравнением (7.110), соответственно для точек ?, = 0 и § = 1 получаем первое и второе дополнительные граничные условия

Для нахождения третьего дополнительного граничного условия продифференцируем уравнение (7.110) по переменной ?, и запишем полученное соотношение для точки Ъ>=

Сравнивая соотношения (7.123) и (7.126), получаем третье дополнительное граничное условие

Используя основные (7.111)—(7.113) и дополнительные (7.124), (7.125), (7.127) граничные условия, можно определить уже шесть коэффициентов ряда (7.115). Подставляя (7.115) во все перечисленные граничные условия, относительно неизвестных коэффициентов Ьк = 0, 5) получаем систему шести алгебраических линейных уравнений. Определяя из решения этой системы коэффициенты Ьк и подставляя их в соотношение (7.115), получаем соотношение

Подставляя (7.128) в (7.117), для определения неизвестной функции q-2.(Po) будем иметь неоднородное обыкновенное дифференциальное уравнение второго порядка

Общее решение уравнения (7.129) разыскивается в виде суммы двух функций:

где ц — частное решение неоднородного уравнения (7.129); ф — общее решение соответствующего однородного уравнения.

Характеристическое уравнение для однородного уравнения будет иметь вид

Последнее уравнение можно представить в виде , { 23,545455 ,

гдеЛ = — матрица собственных значении характеристиче-

' Г1 о!

ского уравнения (7.131); Е = — единичная матрица.

Координаты собственных векторов матрицы Л, отвечающих собственному значению р, удовлетворяют однородной системе уравнений (Л + рХ) = О, которая в данном случае будет иметь вид

371

Решением этой системы уравнений является уравнение вида (7.131), из которого находятся собственные числа р] и р2, имеющие значения pj = = -2,471; Р2 = -22,074. Полученные собственные числа незначительно отличаются от собственных чисел краевой задачи Штурма—Лиувилля, точные значения которых равны pt = 2,467401; р2 = 22,206609.

Следовательно, рассмотренный выше метод решения задачи во второй стадии процесса эквивалентен решению краевой задачи Штурма —Лиувилля.

Общее решение однородного дифференциального уравнения принимает вид

где С и С2 — константы интегрирования.

Частное решение неоднородного уравнения (7.129) принимается в виде г) = Сз, где константа С3 находится из выполнения уравнения (7.129). В данном конкретном случае С3 = 1.

С учетом найденных частного решения неоднородного и общего решения соответствующего однородного уравнения соотношение (7.130) примет вид

Константы интегрирования С и С2 находятся из начальных условий q2(Fox) = 0; q2(Fo{)/dFo = 0.

Формулы для них будут иметь вид

С учетом С и С2 формула для q2(Fo) приводится к виду q2(Fo) = 1 - l,1261exp[-2,4709(Fo - Fo{)] + 0,1261ехр[-22,0745(Fo - Fo,)|. (7.132)

Подставляя (7.132) в (7.128), находим окончательное выражение для решения задачи (7.110)—(7.113) во втором приближении второй стадии процесса:

где Fo, = 0,05 (найдено во втором приближении первой стадии процесса).

Результаты расчетов безразмерных температур по формуле (7.63) в сравнении с точным решением [49] представлены на графиках рис. 7.15. Их анализ позволяет заключить, что полученное здесь решение во всем диапазоне изменения переменной Фурье второй стадии процесса практически совпадает с точным.

Анализируя невязки уравнения (7.110), можно заметить, что во втором приближении второй стадии процесса в диапазоне времени Fo, < Fo < °° это уравнение удовлетворяется практически точно (рис. 7.16, 7.17).

Ввиду достаточно высокой точности решения во втором приближении второй стадии процесса необходимость выполнения третьего приближения имеет лишь теоретическое значение. Для получения решения задачи (7.110)—(7.113) в третьем приближении найдем следующие три дополнительных граничных условия. Для этого продифференцируем по переменной Фурье дополнительные граничные условия (7.124), (7.125), (7.127).

Изменение невязки е уравнения (7.110) по координате ? при Fo = 0,5

Рис. 7.16. Изменение невязки е уравнения (7.110) по координате ? при Fo = 0,5:

------первое приближение;

о — второе приближение

Изменение невязки г уравнения (7.110) для ? = 0,6

Рис. 7.17. Изменение невязки г уравнения (7.110) для ? = 0,6:

------первое приближение;

о — второе приближение

Сравнивая полученные соотношения с соотношениями (7.94), (7.95), применительно к точкам ? = 0 и ?, = 1 будем иметь соотношения

Подставляя (7.115) в основные (7.111)—(7.113) и дополнительные (7.124), (7.125), (7.127), (7.134) граничные условия, относительно неизвестных коэффициентов bk (k = 0, 8) будет иметь систему девяти алгебраических линейных уравнений. Ее решение имеет вид

где v = 1 - q% V = 17/28; v2 = 19/1680; t>3 = 17/4; v4 = 17/24; v5 = 37/4; v6 = = 5/16; v1 = 25/2; г;8 = 7/15; = 27/4; v10 = 13/48; vn = 19/14; v12 = 2/35.

Подставляя найденные значения b/, в (7.115), получаем соотношение

Подставляя (7.135) в (7.117), относительно неизвестной функции q2{Fo) получаем обыкновенное дифференциальное уравнение третьего порядка

Граничные условия для уравнения (7.136) имеют вид

Интегрируя (7.136), получаем уравнение

где Су С С3 — постоянные интегрирования, определяемые из граничных условий (7.137). Формулы для них будут иметь вид

373

С учетом найденных значений С1? С2, С3 формула (7.138) принимает вид

Соотношения (7.135), (7.139) представляют решение задачи (7.110)—

(7.113) в третьем приближении. Результаты расчетов по формуле (7.135) в сравнении с точным решением [49] даны в табл. 7.6.

Таблица 7.6

4

Fo = 0,035

Fo = 0,4

Fo = 1,0

0(? 0,035)

0(4; 0,4)

0(4; 1-0)

По формуле (7.135)

Точное

решение

По формуле (7.135)

Точное

решение

По формуле (7.135)

Точное

решение

0,2

0,5474

0,5503

0,1469

0,1467

0,0334

0,0333

0,4

0,8734

0,8694

0,2796

0,2789

0,0635

0,0634

0,6

0,9835

0,9767

0,3847

0,3839

0,0875

0,0873

0,8

0,9996

0,9975

0,4522

0,4513

0,1029

0,1027

1,0

0,9999

0,9997

0,4755

0,4745

0,1082

0,1079

Отметим, что в третьем приближении произошло существенное уточнение первого и второго собственных чисел в соотношении (7.139). Точное значение третьего собственного числа р3 = 61,6850.

Характерной особенностью полученных здесь аналитических решений является полиномиальная зависимость температуры от координаты ^ в отличие от классических точных аналитических решений, где такая зависимость выражается через тригонометрические функции. Полиномиальная зависимость позволяет получить решение в виде поля изотермических линий. Принцип построения изотерм рассмотрим на примере первого приближения первой и второй стадий процесса. Выражая координату i; как функцию температуры ©(?,, Fo) и времени Fo> соотношения (7.76) и (7.120) можно привести к виду

где Е = ехр[-3(/го - Fot)].

Соотношения (7.140), (7.141) позволяют для любых конкретных 0(^, Fo) — const построить графики зависимости температур от ?, и Fo (графики изотерм, рис. 7.18).

Отметим, что нулевая изотерма 0(^, Fo) = 0 совпадает с графиком движения фронта температурного возмущения по координате в зависимости от времени Fo в первом приближении (см. рис. 7.13). В самом деле, при 0(%, Fo) = 0 выражение (7.140) принимает вид соотношения (?, = V 127b), полностью совпадающего с формулой (7.79), характеризующей перемещение фронта температурного возмущения. Отсюда следует, что в физическом смысле фронт температурного возмущения является аналогом изотермы, движущейся по координате во времени Fo. В данном случае это есть нулевая изотерма — изотерма начального условия. Ввиду положенной в ос-

Распределение изотерм @(?, Fo) = const в координатах Fo(первое приближение, Fo = 0,0833)

Рис. 7.18. Распределение изотерм @(?, Fo) = const в координатах Fo(первое приближение, Fot = 0,0833)

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

Первые производные по времени от соотношений (7.140), (7.141) позволяют определить безразмерные скорости движения изотерм v = d^/dFo по координате ^ в зависимости от времени, а вторые производные — ускорения а = d2Z)/dFo2. Формулы скоростей для первой и второй стадий процесса соответственно будут равны

Графики скоростей движения изотерм, найденных по формулам (7.142),

(7.143), даны на рис. 7.19. Их анализ позволяет заключить, что максимальные скорости изотермы имеются вблизи точки ? = 0. По мере удаления от этой точки скорости существенно уменьшаются, достигая некоторого минимума. Затем при приближении к точке ? = 1 скорости изотерм вновь значительно возрастают. Графики ускорений изотерм по форме (качественно) практически соответствуют графикам скоростей и отличаются от них лишь количественно.

Ввиду невысокой точности первого приближения (и особенно во второй стадии процесса) изотермы, определяемые по формулам (7.140), (7.141), имеют небольшой излом при Fo = Fo = 0,0833, т.е. в точке сопряжения решений для первой и второй стадий процесса (см. рис. 7.18). В связи с чем на графиках рис. 7.19 в этой точке имеет место некоторый скачок в эпюрах скоростей, который уже во втором приближении практически не наблюдается, так же как и излом в изотермах (рис. 7.20, 7.21).

Графики скоростей изотерм (первое приближение)

Рис. 7.19. Графики скоростей изотерм (первое приближение)

Распределение изотерм Рис. 7.21. Изменение скоростей изотерм во втором приближении во втором приближении

Рис. 7.20. Распределение изотерм Рис. 7.21. Изменение скоростей изотерм во втором приближении во втором приближении

Чтобы построить изотермические линии в координатах Fo применительно к последующим приближениям, для каждых конкретных ©(%, Fo) и Fo относительно % необходимо решать алгебраическое уравнение. Ввиду того что каждому ?, и Fo согласно полученным аналитическим решениям соответствует лишь одно значение температуры 0(?, Fo), алгебраическое уравнение имеет лишь один корень, удовлетворяющий соответствующим решениям вида (7.90), (7.102), (7.105), (7.108), (7.109), (7.133), (7.135).

Объяснение распределению скоростей, показанному на рис. 7.20, 7.21, можно получить из анализа формул (7.142), (7.143). И, в частности, из формулы (7.142) следует, что с уменьшением числа Фурье (Fo -* 0) скорости изотерм неограниченно возрастают. Анализ формулы (7.143) позволяет заключить, что с приближением числа Фурье к значению, при котором соответствующая изотерма достигает координаты ?, = 1 (см. рис. 7.18), скорость изотермы также неограниченно возрастает. Исходя из формулы (7.143), неограниченное возрастание скорости может быть лишь в случае, когда знаменатель этой формулы приближается к нулю. В самом деле, определяя, например, для 0 = 0,3 величину числа Фурье, при котором знаменатель в формуле (7.143) обращается в нуль, получаем Fo = 0,20213. Расчеты показывают, что именно при этом числе Фурье изотерма 0 = 0,3 достигает центра пластины (см. рис. 7.18). Аналогичная ситуация имеет место и для любых других изотерм.

Найдем тепловой поток, приходящийся на единицу площади ограничивающей поверхности пластины:

Из последнего соотношения следует, что величина теплового потока прямо пропорциональна коэффициенту теплопроводности тела, разности между температурой стенки и начальной температурой и обратно пропорциональна Vx. Следовательно, в начальный момент времени тепловой поток бесконечно велик. Такая величина теплового потока приводит к бесконечно большим скоростям движения изотерм. С увеличением времени тепловой поток уменьшается, приводя к соответствующему уменьшению скоростей изотерм. Кроме того, при малых значениях времени (т —> 0) величина прогретого слоя совершенно незначительна (q(Fo) —? 0) и его термическое сопротивление практически отсутствует (t](Fo)/X —> 0), что является причиной бесконечно большой величины теплового потока.

 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ     След >