Математическая модель сушки коры деревьев при высокоинтенсивном нагреве Текст научной статьи по специальности «Физика»

4. Tulloch, M. Windows 7 - Resource Kit / Mitch Tul- ton 98052-6399: Microsoft Press, 2010. - P. 90 - 91. loch, Tony Northrup, Jerry Honneycutt-Redmont. - Washing-

УДК 519.63

Н.Н. Синицын, З.К. Кабаков, Д.А. Домрачев

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ СУШКИ КОРЫ ДЕРЕВЬЕВ ПРИ ВЫСОКОИНТЕНСИВНОМ НАГРЕВЕ

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

Математическая модель сушки коры деревьев, тестирование, задача Стефана, погрешность решения.

The article presents the mathematical description of the process of the bark drying at high-intensity warming up and the way of testing of a numerical model of bark drying using the exact solution of the Stefan problem. The results of influence of adjusting parameters of a numerical algorithm on the modeling error are presented in the paper.

Mathematical model of bark drying, testing, Stephane's task, decision error.

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

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

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

■ начальное условие:

ЭТ Э

ЭГ

сЭф (T) P(T) — = — (1(T) — ), dt Эх Эх

интегрируемое в области: 0 < х < S, 0 < t < tk ;

(1)

Tl t=0 T

- граничное условие:

ЭТ

при х = 0 -1(Т)— = a(T - T );

Эх

при

х = S

1(T ) — = 0, Эх

(2)

(3)

(4)

где р - плотность материала; а - коэффициент теплоотдачи; - половина толщины пластины; 71 -

температура среды; 7 - начальная температура материала; 1 - коэффициент теплопроводности.

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

' С (Т), Т > Тл ; с(Т с ) '¥+ с(Т л ) • (1 -y) + С2(Т), Т < Тс,

gjJ^ Т < T < T ■

DT ' с л'

коэффициент теплопроводности и плотность определяют по формулам:

1 =

Ii, Т < Тс;

1 •y+12(1 -У), Тс < T < T ;

l2, Т > Тл,

6эф =

Р =

Т > Тл;

Рl,

Р1 у + Р2(1 -У), ТС < Т < Тл:

х1 = (, - 0,5) Ах , для дискретных моментов времени

Р2

Т < Т

где Тл = Тф +19,5, Тс = Тф - 37 - фиктивные температуры начала и окончания фазового перехода воды; с (Т), - теплоемкость материала; с1 и с2 - теплоемкость сухого и влажного слоев материала; р1 и р2 - плотность сухого и влажного слоев материала; 11 и 12 - коэффициенты теплопроводности сухого и влажного слоев материала; g - доля влаги в элементарном объеме материала; - половина толщины слоя материала; Ь - удельная теплота фазового перехода влаги; у - доля влажного материала.

Величина у определяется по формуле:

1, Т < Т„;

у = ^

Т л - Т Т л - Т„

Т < Т < Т •

0, Т > Тл.

На рис. 1 показана схема расчетной области.

Рис. 1. Схема расчетной обл асти Расчетная область: 1 - сухая зона; 2 - влажная зона; 3 - двухфазная зона; А1 - ширина двухфазной зоны; ес, е, ел - координаты границ начала двухфазной зоны, фазового перехода, окончания двухфазной зоны соответствующих температур Тс, Тф,

Тл и Тп - температуры поверхности.

Система уравнений (1) - (4) в общем случае может быть решена только численным методом. При использовании метода конечных разностей (МКР) значение температур определяют в узлах расчетной области, координаты которых находят по формуле:

хп = Ах • п , где , = 0, N +1, N - количество узлов внутри расчетной области, 0 и N +1 - номера фиктивных узлов, находящихся за пределами области на расстоянии Ах /2; Ах = 5 / N - расстояние между

узлами; п = 0, [хк / Ах] - моменты времени (п = 0 -начальный момент времени); Ах - расчетный шаг по времени. Для краткости температуры Т(х1, х") -

обозначают Т".

При использовании явной схемы аппроксимации производных по координате температуру в следующий момент времени п + 1 в N внутренних узлах определяют по формуле:

Т П+1 _ Тп +

Ах

с(Тп) • Р(Тп) •Ах2

1+1 • (Т+1 - Т) -1,. -! • (Т - Т-,)

грп . грп грп грп

где , = , 1 , =1 -'-, 1 , =1 т--'-.

1+1 2 1-1 2 2 2

Температуру в начальный момент времени задают по формуле:

Т1 = Т0 для , = 0, N +1.

Температуру в фиктивных узлах: , = 0 и N +1 в момент времени п + 1 определяют по формулам:

Т= 1 0

(1 -с) Т + 2 •%• Т, 1+х

% = -

а-Ах 21 ,

Т = Т

Т N+1 N'

Расположение границы перехода воды в пар определяют в поле температуры по температуре фазового перехода влаги в цикле по , = 2..^ из условия:

если Т-1 > 7ф > Т;

то е = Ах(, - 3/2) + Ах

Т - Т Т - Т

Численное решение при явной схеме аппроксимации является условно устойчивым. В этом случае расчетный шаг определяем по формуле: Ах = Ах2 / (ку • а), где ку > 2.

Погрешность численного решения будет зависеть от настроечных параметров алгоритма N ку и АТ .

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

Для выбора этих параметров выполним тестирование численного решения задачи Стефана путем сравнения с точным решением этой задачи [1], кото-

п

х

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

Т,(х,х) = Тп + (Тф -Тп)

Г (^=)

2у] а1х

ф Т п> р

ег/ ) 2 а1

(5)

- во влажной зоне:

ег/с( .-)

Т2 (х, х) = Т0 - (Т0 - Тф)--^,

ег/с( ) 2 а2

(6)

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

е

(7)

где 0 - корень трансцендентного уравнения (8): ^ 12 • (70- ^

-х

11(Тф - Т0) р2 12 • (Т0 - Тф)

4а • ег/ (т"г=)

^ а

р ехр(-^-) -Р ч 4а,

1 л/а2 • ег/с(^=)

Р

х ехр(—-) = р2 • я • I Р—.

(8)

В формулах (5) - (8) использованы следующие обозначения: 11 и 12 - теплопроводности сухой и влажной зон материала, а1 = 1, / (с,. • р,.) - температуропроводность (I = 1 - сухая зона, . = 2 - влажная зона), с1, р,. - теплоемкость и плотность, Тп - температура поверхности, Т0 - начальная температура ма-!_

о

териала,

ег/ (х) = ^ | е"

ёX, ег/с(х) = 1 - ег/(х) =

О ^

= А } ёX .

х

Выполним тестирование для пластины толщиной 55 = 0,0125 м, нагреваемой в симметричных условиях. Исходные данные для моделирования и расчета по формулам (5) - (8) приведены в таблице.

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

Из рис. 2 видно, что при N = 140 численное решение практически совпадает с точным.

4а

2

Таблица

Исходные данные

№ п/п Величина, размерность Значение величины

В модели В точном значении

1 Половина толщины, 5 м 0,0125 ¥

2 Начальная температура, °С 180 180

3 Температура поверхности, °С 20

4 Температура среды, °С 20

5 Плотность влажной зоны, кг/м3 821 821

6 Плотность сухой зоны, кг/м3 718 718

7 Коэффициент теплопроводности влажной зоны, Вт/(м- К) 0,4629 0,4629

8 Коэффициент теплопроводности сухой зоны, Вт/(м- К) 0,4629 0,4629

9 Теплоемкость влажной зоны, Дж(кг-К) 3498 3498

10 Теплоемкость сухой зоны, Дж(кг-К) 2720 2720

11 Температуры фазового перехода, °С 100 100

12 Удельная теплота фазового перехода, Дж/кг 2256800 2256800

13 Коэффициент теплоотдачи, Вт(м2 К) 1010

14 Конечное время процесса, с 245,36 245,36

15 Фиктивный интервал ДТ, °С 56,5

16 Доля влаги, кг/кг сухого зоны 0,25 0,25

т , °С 160 ■ 140 ■ 120 100 30 ■ 60 ■ 40 20 0 ■

\

\ •

\

\

X

--

А

\ , 1

10

X • 103

Рис. 2. Распределение температуры по толщине пластины 1 - решение по модели; 1' - точное решение; модель при N = 40; к = 3;

Т л = 119,5 °С ; Тс

63 °С; ДТ = 119,5 - 63 = 56,5 °С; х = 245,36 °С

м

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

о 11

§ = = А -VI (е.-еП )2-100%,

где п = хк / Ах , е - среднеарифметическое значение координаты фазового перехода воды; еп - результат моделирования координаты фазового перехода в момент времени п, еп - точное решение в момент времени п.

Результаты моделирования приведены на рис. 3.

Как видно из рис. 3а, с увеличением N более 640 узлов погрешность асимптотически стремится к нулю. При увеличении ку погрешность уменьшается.

При ку =15 погрешность начинает увеличиваться.

Учитывая, что при ку = 3 погрешность меньше 1 % -

в дальнейшем исследовании принимаем ку = 3.

Исследование влияния фиктивного интервала ДТ показало, что имеется несимметричный интервал «размазывания» теплоемкости в окрестности температуры фазового перехода, при котором численное решение по динамике координаты перехода имеет наименьшее значение среднеквадратичной погрешности.

Следует отметить, что такое большое количество узлов, обеспечивающее погрешность менее 1 %, обусловлено «жестким» граничным условием I рода. При исследовании граничных условий III рода следует ожидать существенного уменьшения допустимого количества узлов.

8,%

1.3

1Д ■

0,9 ■

0,7 "

0,5 ■

0,3 ■

0,1 ■

8, % 1,4 1,2 1

0,6

120 юз гас зео адо 520 воо

N

а)

\

\ \

\ \ \

ч

0,4

б)

Рис. 3. Зависимость средней квадратичной погрешности 8 от количества узлов N и параметра ку при параметре ДТ в интервале температур 63 ^ 119,5 °С. а - параметр ДТ в интервале 63 ^ 119,5 °С, ку = 2; б - параметр ДТ в интервале 63 ^ 119,5 °С, N = 40.

п=1

к

у

Увеличение количества узлов и соответствующее уменьшение шага по времени, согласно условию устойчивости, влияет на уменьшение погрешности более эффективно, чем только уменьшение расчетного шага по времени при увеличении ку .

В результате тестирования установлено, что для уменьшения средней относительной погрешности до 1 % необходимо взять количество узлов сетки не бо-

лее 40. Следует отметить, что для менее «жесткого» граничного условия, например, III рода, это ограничение может измениться до существенно меньших N.

Литература

1. Лыков, А.В. Теория теплопроводности / А.В. Лыков. - М., 1967.