Scientific journal
Modern high technologies
ISSN 1812-7320
"Перечень" ВАК
ИФ РИНЦ = 0,940

MATHEMATICAL MODELING OF THERMAL PROCESSES IN THE TUBULAR ELECTRIC HEATER OF THE BIOREACTOR

Ozhogova E.V. 1 Lubentsova E.V. 1 Lubentsov V.F. 1
1 Kuban State Technological University
Taking into account the dynamics of heat flows that create a perturbation to the environment temperature through the surface of the heater wall, a mathematical model of the heater used to regulate the temperature of the environment of the anaerobic fermentation process in a bioreactor is formalized. The model is based on the equations of heat transfer through the wall using the classical method of solving the thermo physics problem – the Fourier method under boundary conditions for the case of a convective heat transfer mechanism. At the same time, the Bio criteria characterizing the unevenness of the temperature field and the Fourier criterion characterizing the unsteadiness of the processes are taken into account in the equations. The calculated ranges of values of the massiveness parameters are given depending on the coefficients of heat transfer by convection in stationary mode from the working environment and from the environment for different values of the heater wall thickness and the coefficient of thermal conductivity of the isotropic wall material. The practical value of the work consists in the fact that, using the obtained ranges of values of the massiveness parameters, the values of the time constant of the model are obtained with a convective heat transfer mechanism that characterizes inertia when temperature is regulated using a tubular heater. The calculated values allow, knowing the properties of the wall material, to estimate the rational thickness of the heater wall, providing minimal inertia at a fixed heating intensity of the environment. In the case of a fixed wall thickness, you can choose the type of heater made of a material with certain (required) thermo physical properties. The model obtained in the work is used in the synthesis and analysis of the temperature control system of a tubular electric heater in a bioreactor.
mathematical model
heater
heat transfer
Bio criteria
time constant

Важная роль в современном развитии энергетической отрасли России принадлежит биоэнергетике, основой которой является использование различных биоматериалов в качестве топлива. Обусловлено это многими факторами, среди которых истощение ископаемых энергоресурсов, а также возможность с помощью биоэнергетики решить экологические проблемы. Биогазовые установки (БГУ) обеспечивают переработку органических отходов (стоков животноводческих производств и растениеводства) и осадков сточных вод в биогаз (горючий газ). Попутно с биогазом на БГУ получают высокоэффективное дорогостоящее жидкое органическое удобрение [1].

Для достижения максимальной эффективности образования биогаза и получения органических удобрений в БГУ должен поддерживаться оптимальный для данной установки температурный режим – важнейший фактор процесса сбраживания [2]. Для обеспечения оптимального температурного режима, а также для лучшего обеззараживания сырья, используются два метода подогрева: прямой подогрев в форме пара или смешивающейся с сырьем горячей воды и непрямой подогрев через теплообменные устройства с различными нагревательными устройствами [3]. Добавление горячей воды повышает влажность субстрата и поэтому может использоваться только там, где это допустимо. Внешний подогрев с помощью теплообменника с теплопроводящими элементами на поверхности стен реактора биогазовой установки менее эффективен из-за потерь тепла с поверхности стен. Непрямой подогрев с размещением теплообменных устройств в зоне действия перемешивающего устройства помогает избежать осаждения твердых частиц на их поверхности. В современных мощных электронагревательных приборах в качестве основы их нагревательных устройств используются трубчатые электрические нагреватели – ТЭНы.

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

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

Материалы и методы исследования

Тепловые потоки свойственны теплообменным агрегатам как периодического, так и непрерывного действия. Функционирование теплообменного агрегата метантенка БГУ проходит в нестационарном тепловом режиме, который характеризуется предварительным нагревом сырья заданного объема и в течение определенного времени перед загрузкой в биореактор и режимом стабилизации температуры среды в биореакторе с высокой точностью, особенно для термофильного режима сбраживания. Решение такой задачи требует получения модели по каналу изменения температуры с учетом динамики тепловых потоков, создающих возмущение на температуру среды через поверхность стенки нагревателя. Возмущение от теплового потока нагревателя влияет на температуру среды в аппарате через стенку нагревателя. Теплопередающая поверхность нагревателя изготавливается из различных материалов, теплофизические свойства которых также изменяются при изменении температуры.

В работе [4] получены основные уравнения переноса тепла через стенку с применением классического метода решения задачи теплофизики – метода Фурье при граничных условиях для случая конвективного механизма теплоотдачи с помощью учета в уравнениях констант – критериев Био, характеризующих неравномерность температурного поля, и критерия Фурье, характеризующего нестационарность процессов [5].

В данной работе с учетом сложных теплофизических эффектов и полученных в [4] основных уравнений теплопереноса через стенку в динамическом режиме формализована математическая модель процесса нагрева среды по каналу возмущения от теплового потока нагревателя (трубчатого электронагревателя, например ТЭНа, теплообменника и др.) для целей моделирования системы управления тепловым агрегатом – метантенком БГУ. Оценка инерционности теплового объекта проведена путем анализа передаточной функции процесса при ступенчатом изменении теплового потока в стенке нагревателя для практического использования этой информации при математическом моделировании динамических процессов в трубчатом электронагревателе.

При конвективном механизме теплоотдачи коэффициентами модели являются коэффициент теплопроводности изотропного материала стенки, полная толщина стенки, а также коэффициенты теплоотдачи конвекцией в стационарном режиме со стороны рабочей среды (α1) и со стороны окружающего пространства (α2). Учет влияния излучения производится суммированием коэффициентов теплоотдачи – конвективного и лучистого. При конвективном механизме теплоотдачи числа массивности стенки æ1 и æ2 , входящие как константы в динамические уравнения, являются числами Био [4]: æ1 = Вi1 = α1 ∙ Δх/λ, æ2 = Вi2 = α2 ∙ Δх/λ, где ∆x – толщина стенки; λ – коэффициент теплопроводности изотропного материала стенки; α1, α2 – коэффициенты теплоотдачи конвекцией в стационарном режиме со стороны рабочей среды и со стороны окружающей среды. Теплопередающие поверхности нагревателей изготавливаются из различных материалов, теплофизические свойства которых также изменяются с изменением температуры. Если представить изменения коэффициентов теплоотдачи от естественной конвекции до конвекции при фазовых переходах в интервале α1 = 10…2∙103 Вт/м∙К, а изменение теплопроводности материалов нагревателя – в интервале λ = 0,01… 400 Вт/м∙К, то область значений числа æ1 при изменении ее толщины Δх = 0,005…0,1 м составляет æ1 = 5 … 0,0025 (при α1 = 10 Вт/м∙К) и æ1 = 1000 … 0,5 (при α1 = 2000 Вт/м∙К). При тех же данных и при изменении коэффициента теплоотдачи α2 в 10 … 400 Вт/м∙К область значений æ2 составляет æ2 = 5 … 0,0025 (при α2 = 10 Вт/м∙К) и æ2 = 200 … 0,1 (при α2 = 400 Вт/м∙К).

Результаты исследования и их обсуждение

В табл. 1, 2 приведены расчетные значения параметров æ1 и æ2 в зависимости от α1 и α2 для различных значений ∆x и λ. При анализе полученных значений следует учитывать, что значения числа Био меньше 0,1 означают, что теплопроводность внутри стенки намного быстрее, чем конвекция тепла вдали от его поверхности, а градиенты температуры незначительны внутри него. Это может указывать на применимость (или неприменимость) в решении данной задачи методов анализа сосредоточенных систем) [6]. Обычно этот метод анализа приводит к простому экспоненциальному нагреванию, поскольку количество тепловой энергии (в широком смысле количество «тепла») в теле прямо пропорционально его температуре, которая, в свою очередь, определяет скорость передачи тепла в него или из него. Это приводит к простому дифференциальному уравнению первого порядка, описывающему теплопередачу в этих системах, что подтверждено полученной моделью в данной работе.

Таблица 1

Значение параметра æ1 в зависимости от α1 при Δх = 0,005; 0,1; λ = 0,01… 400

№ 1

α1

λ = 0,01

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

5

25

50

100

200

300

400

500

1000

 

Δх = 0,1

æ1

100

500

1000

2000

4000

6000

8000

10000

20000

№ 2

α1

λ = 0,1

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,5

2,5

5

10

20

30

40

50

100

 

Δх = 0,1

æ1

10

50

100

200

400

600

800

1000

2000

№ 3

α1

λ = 10

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,005

0,025

0,05

0,1

0,2

0,3

0,4

0,5

1

 

Δх = 0,1

æ1

0,1

0,5

1

2

4

6

8

10

20

№ 4

α1

λ = 20

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,0025

0,0125

0,025

0,05

0,1

0,15

0,2

0,25

0,5

 

Δх = 0,1

æ1

0,05

0,25

0,5

1

2

3

4

5

10

№ 5

α1

λ = 40

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,0013

0,006

0,013

0,025

0,05

0,075

0,1

0,125

0,25

 

Δх = 0,1

æ1

0,025

0,125

0,25

0,5

1

1,5

2

2,5

5

№ 6

α1

λ = 50

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,001

0,005

0,01

0,02

0,04

0,06

0,08

0,1

0,2

 

Δх = 0,1

æ1

0,02

0,1

0,2

0,4

0,8

1,2

1,6

2

4

№ 7

α1

λ = 100

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,0005

0,0025

0,005

0,01

0,02

0,03

0,04

0,05

0,1

 

Δх = 0,1

æ1

0,01

0,05

0,1

0,2

0,4

0,6

0,8

1

2

№ 8

α1

λ = 150

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,000333

0,0016

0,0033

0,006

0,013

0,02

0,026

0,033

0,067

 

Δх = 0,1

æ1

0,006667

0,0333

0,0666

0,133

0,267

0,4

0,533

0,667

1,333

№ 9

α1

λ = 200

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,00025

0,00125

0,0025

0,005

0,01

0,015

0,02

0,025

0,05

 

Δх = 0,1

æ1

0,005

0,025

0,05

0,1

0,2

0,3

0,4

0,5

1

№ 10

α1

λ = 250

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,0002

0,001

0,002

0,004

0,008

0,012

0,016

0,02

0,04

 

Δх = 0,1

æ1

0,004

0,02

0,04

0,08

0,16

0,24

0,32

0,4

0,8

№ 11

α1

λ = 300

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,00017

0,0008

0,0016

0,003

0,007

0,01

0,0133

0,017

0,033

 

Δх = 0,1

æ1

0,00333

0,0166

0,0333

0,066

0,133

0,2

0,266

0,333

0,667

№ 12

α1

λ = 400

10

50

100

200

400

600

800

1000

2000

 

Δх = 0,005

æ1

0,00013

0,000625

0,00125

0,0025

0,01

0,01

0,01

0,01

0,03

 

Δх = 0,1

æ1

0,0025

0,0125

0,025

0,05

0,1

0,15

0,2

0,25

0,5

Таблица 2

Значение параметра æ2 в зависимости от α2 при Δх = 0,005; 0,1; λ = 0,01 … 400

№ 1

α2

λ = 0,01

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

5

10

15

20

30

40

50

100

200

 

Δх = 0,1

æ2

100

200

300

400

600

800

1000

2000

4000

№ 2

α2

λ = 0,1

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,5

1

1,5

2

3

4

5

10

20

 

Δх = 0,1

æ2

10

20

30

40

60

80

100

200

400

№ 3

α2

λ = 10

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,005

0,01

0,015

0,02

0,03

0,04

0,05

0,1

0,2

 

Δх = 0,1

æ2

0,1

0,2

0,3

0,4

0,6

0,8

1

2

4

№ 4

α2

λ = 20

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,0025

0,005

0,008

0,01

0,015

0,02

0,025

0,05

0,1

 

Δх = 0,1

æ2

0,05

0,1

0,15

0,2

0,3

0,4

0,5

1

2

№ 5

α2

λ = 40

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,001

0,003

0

0,01

0,01

0,01

0,01

0,03

0,05

 

Δх = 0,1

æ2

0,025

0,05

0,08

0,1

0,15

0,2

0,25

0,5

1

№ 6

α2

λ = 50

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,001

0,002

0,003

0,004

0,006

0,008

0,01

0,02

0,04

 

Δх = 0,1

æ2

0,02

0,04

0,06

0,08

0,12

0,16

0,2

0,4

0,8

№ 7

α2

λ = 100

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,0005

0,001

0,002

0,002

0,003

0,004

0,005

0,01

0,02

 

Δх = 0,1

æ2

0,01

0,02

0,03

0,04

0,06

0,08

0,1

0,2

0,4

№ 8

α2

λ = 150

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,00033

0,000667

0,001

0,001

0,002

0,003

0,003

0,007

0,013

 

Δх = 0,1

æ2

0,00667

0,013333

0,02

0,027

0,04

0,053

0,067

0,133

0,267

№ 9

α2

λ = 200

10

20

30

40

60

80

100

200

400

 

Δх =0,005

æ2

0,00025

0,0005

0,00075

0,001

0,0015

0,002

0,0025

0,005

0,01

 

Δх = 0,1

æ2

0,005

0,01

0,015

0,02

0,03

0,04

0,05

0,1

0,2

№ 10

α2

λ = 250

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,0002

0,0004

0,0006

0,0008

0,0012

0,0016

0,002

0,004

0,008

 

Δх = 0,1

æ2

0,004

0,008

0,012

0,016

0,024

0,032

0,04

0,08

0,16

№ 11

α2

λ = 300

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,00017

0,00033

0,0005

0,00067

0,001

0,001

0,0017

0,003

0,0067

 

Δх = 0,1

æ2

0,00333

0,00667

0,01

0,013

0,02

0,027

0,033

0,067

0,133

№ 12

α2

λ = 400

10

20

30

40

60

80

100

200

400

 

Δх = 0,005

æ2

0,00013

0,00025

0,00038

0,0005

0,0008

0,001

0,0013

0,0025

0,005

 

Δх = 0,1

æ2

0,0025

0,005

0,0075

0,01

0,015

0,02

0,025

0,05

0,1

В общем случае влияние изменяющейся температуры стенки на тепловой поток описано в работе [4] передаточной функцией, содержащей сумму бесконечного ряда, следующего вида:

W*(p) = k / (1 + a1 р + a2 р2 + … ) = k / (1+aj ∙ рj), j =1, 2 … . (1)

где k = æ2 / (æ1 + æ2 + æ1 ∙ æ2) – коэффициент передачи;

missing image file, j =1, 2 … .

Ограничившись в выражении (1) первыми членами (приняв j = 1), как правило, наиболее существенными в сумме, получим модель влияния теплового потока на температуру стенки в виде передаточной функции апериодического звена первого порядка:

W(p) = k / (1+ Т∙р), (2)

где T = α1 ∙ Ts – постоянная времени, определяемая соотношениями

missing image file;

Ts = α1 ∙ ρ ∙ c ∙ (Δх)2 / λ; λ – коэффициент теплопроводности; ρ – плотность; с – удельная массовая теплоемкость.

Это приближение мало отличается от параметров переходных характеристик, рассчитанных с помощью сумм бесконечных рядов. Передаточная функция (2) соответствует процессу переноса тепла в динамическом режиме в области изображений по Лапласу относительно безразмерного времени, т.е. числа Фурье Fo = λ∙τ/ρ∙c (∆x)2. Для получения передаточной функции в области изображений относительно размерного времени τ был заменен аргумент-оператор р на Ts∙p, где Ts – постоянная времени, определяемая соотношением Ts=ρ∙c∙(∆x)2/λ. Эта замена без сложных преобразований объясняется тем, что в определении изображений по Лапласу функции комплексного переменного F(p)=∫f(τ)∙e-pτ dτ независимые постоянные ∆x и λ/ρ∙c выносятся за знак интеграла [7].

Физические характеристики λ, ρ, c функционально связаны с коэффициентом температуропроводности, характеризующим темп нагрева тела, соотношением а = λ /ρ∙c [8]. Коэффициент температуропроводности для стальной стенки можно охарактеризовать значением а ≈ 10-5 м2/с. В нашем случае при расчете Тs используется величина, обратная коэффициенту а (1/а=ρ∙c/λ), являющаяся тепловой инерцией. С учетом значений а и тепловой инерции, равной 105, постоянная времени равна Т = а1∙а∙(∆x)2 = а1 ∙(∆x)2∙105.

Для адаптации полученной модели в качестве аппроксимирующей функции инерционности объекта, например, ТЭНа при ступенчатом изменении теплового потока необходимо рассмотреть конкретные значения постоянной времени, характеризующей динамику изменения температуры в стенке нагревателя при изменении условий теплопереноса. Так, для условий стальной стенки и различных значений коэффициента теплопроводности и коэффициентов æ1 = 0,01 и æ2 в интервале 0,01 … 10,0; æ2 = 0,1 и æ1 в интервале 0,01 … 10,0 при изменении Δх от 0,005 до 0,1 м с использованием рассчитанных в табл. 3 с учетом тепловой инерции, равной 105, значений параметра а1 получены значения постоянной времени Т , приведенные в табл. 4.

Таблица 3

Значения параметра а1 в функции чисел массивности æ1, æ2

 

æ1

0,01

0,05

0,1

0,25

0,5

0,75

1

2,5

5

7,5

10

æ2

0,01

50,25

17,03

9,51

4,31

2,44

1,79

1,48

0,89

0,69

0,63

0,59

0,05

17,03

10,25

6,94

3,60

2,22

1,68

1,39

0,86

0,67

0,61

0,58

0,10

9,51

6,94

5,25

3,14

2,01

1,55

1,31

0,82

0,65

0,59

0,56

0,25

4,31

3,69

3,14

2,24

1,59

1,29

1,11

0,73

0,59

0,54

0,51

0,50

2,44

2,22

2,01

1,59

1,23

1,04

0,92

0,64

0,52

0,48

0,46

0,75

1,79

1,68

1,55

1,29

1,04

0,89

0,80

0,57

0,47

0,44

0,42

1,00

1,48

1,39

1,31

1,11

0,92

0,80

0,72

0,53

0,44

0,41

0,39

2,50

0,89

0,86

0,82

0,73

0,64

0,57

0,53

0,40

0,34

0,32

0,30

5,00

0,69

0,67

0,65

0,59

0,52

0,47

0,44

0,34

0,29

0,27

0,26

7,50

0,63

0,61

0,59

0,54

0,48

0,44

0,41

0,32

0,27

0,25

0,24

10,0

0,59

0,58

0,56

0,51

0,46

0,42

0,39

0,30

0,26

0,24

0,23

Таблица 4

Значения параметра Tq в функции толщины стенки ∆x, а1 при æ1 и æ2, равных 1; 2,5; 5; 7,5; 10,0

 

а1

50,25

10,25

5,25

2,24

1,23

0,89

0,72

0,40

0,29

0,25

0,23

∆x

0,005

125,6

25,6

13,1

5,6

3,08

2,23

1,8

1

0,73

0,63

0,58

0,01

502,5

102,5

52,5

22,4

12,3

8,9

7,2

4

2,9

2,5

2,3

0,05

12562,5

2562,5

1312,5

560

307,5

222,5

180

100

72,5

62,5

57,5

0,1

50250

10250

5250

2240

1230

890

720

400

290

250

230

missing image file

Рис. 1. Зависимость постоянной времени нагревателя в функции толщины стенки Δх при æ1 = æ2 =1,0 … 10

missing image file

Рис. 2. Зависимость постоянной времени нагревателя в функции параметра а1, характеризующего условия теплоотдачи при конвективном механизме

Используя предварительно рассчитанные области значений чисел æ1 и æ2, в табл. 4 приведены расчетные значения постоянной времени нагревателя при конвективном механизме теплоотдачи.

На рис. 1 представлены зависимости изменения постоянной времени Т в функции характерного линейного размера стенки Δх при значениях чисел массивности стенки æ1 и æ2 = 1,0 … 10,0, определяющих значения параметра а1 от 0,22 до 0,72 при конвективном механизме теплоотдачи. Заметно сильное влияние на Т больших значений Δх. При фиксированной толщине стенки Δх постоянная времени Т, как видно из рис. 2, линейно связана с изменениями параметра а1, характеризующего условия теплоотдачи при конвективном механизме.

Заключение

Изменение плотности теплового потока со стороны рабочей или внешней среды следует рассматривать как возмущение теплового потока на объект управления, например на трубчатый теплоэлектронагреватель. Со стороны внешней среды могут происходить изменения, связанные, например, с изменением скорости движения биомассы в биореакторе в результате перемешивания, отражающиеся на величине теплового потока. В этих условиях влияния тепловых потоков со стороны рабочей и внешней среды отражаются на изменении температуры стенки. Используя математическую модель теплопереноса через стенку нагревателя, можно получить значения постоянной времени звена, характеризующего инерционность объекта по каналу регулирования температуры в биореакторе. Зная свойства материала стенки, можно оценить рациональную толщину Δх, обеспечивающую минимальное запаздывание при фиксированной интенсивности нагрева среды. В случае фиксированной Δх можно выбрать тип нагревателя, изготовленного из материала, обладающего определенными (требуемыми) теплофизическими свойствами и предельной температурой.

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