АЛГОРИТМ И ЧИСЛЕННАЯ СХЕМА МОДЕЛИРОВАНИЯ НЕСТАЦИОНАРНЫХ ПРОЦЕССОВ ТЕПЛОВОЙ ОБРАБОТКИ БЕТОННЫХ ИЗДЕЛИЙ Текст научной статьи по специальности «Строительство и архитектура»

УДК 004.94:666.97.035

АЛГОРИТМ И ЧИСЛЕННАЯ СХЕМА МОДЕЛИРОВАНИЯ НЕ СТАЦИОНАРНЫХ ПРОЦЕССОВ ТЕПЛОВОЙ ОБРАБОТКИ БЕТОННЫХ ИЗДЕЛИЙ

А.М. НИЯКОВСКИЙ (Полоцкий государственный университет);

Ю.В. ЯЦКЕВИЧ

(Белорусский государственный университет информатики и радиоэлектроники);

А.Н. ЧИЧКО

(Белорусский национальный технический университет)

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

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

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

Подход, основанный на эмпирических зависимостях. Данный эмпирический подход [1-3], сформировался в рамках науки о бетоне и строительного материаловедения и заключается в экспериментальном подборе режимов ТВО и проектировании оптимальных составов бетонных смесей с последующим обобщением полученных простых факторных зависимостей в виде графиков, таблиц и аппроксимирующих уравнений. Метод оперирует средними по объему обрабатываемого изделия теплофизическими характеристиками материалов, температурой и степенью гидратации и не принимает в расчет начальные, граничные и краевые условия в пределах теплотехнологической установки для проведения ТВО.

Комплексный поход, базирующийся на методах системного анализа широкого спектра технологических параметров. Данный подход использует метод многофакторного проектирования режимов ТВО [4-7] и заключается в получении «математической модели» [4], отражающей влияние различных технологических факторов в их взаимосвязи на требуемые целевые характеристики обрабатываемого изделия, такие как прочность и степень гидратации. Зависимость этих целевых характеристик от отдельных факторов задана уравнениями, полученными экспериментальным путем и обработанными с учетом имеющихся материаловедческих представлений о процессе твердения бетона. В системе уравнений, определяющих такую математическую модель, могут учитываться некоторые теплотехнические параметры, например, теплообмен с окружающей средой [4, 6]. Однако такая модель по-прежнему оперирует интегральными характеристиками и не учитывает особенности нестационарного процесса нагрева и охлаждения бетонного изделия в теплотехнологической установке ускоренной гидратации, поскольку не использует аналитический аппарат математической физики.

Подход, основанный на аналитических методах. Аналитический подход заключается в использовании фундаментальных уравнений математической физики, адаптированных к задачам твердения бетона в условиях подвода теплоты, описывающих нестационарную теплопроводность при наличии внутренних тепловыделений [8-11]. Использование данного метода ограничено трудностями отыскания аналитических решений для тел сложной геометрии при трехмерной постановке задачи в условиях переменных теплофизических свойств нагреваемых материалов и оборудования [9].

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

с учетом всех конструкционных, технологических и теплотехнических факторов, определяющих процессы, протекающие в теплотехнологических установках ускоренной гидратации бетона [11-17].

Целью исследования является разработка основанных на нестационарном трехмерном уравнении теплопроводности, учитывающем источник внутренних тепловыделений, новых алгоритмов компьютерного моделирования энергосберегающих режимов ТВО бетона, позволяющих выполнять численные расчеты процессов нагрева и гидратации для бетонных изделий произвольной 3Б-геометрии, обладающих сложной композитной структурой.

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

Разработанная теплофизическая модель процесса включает систему граничных и начальных условий, различающихся в зависимости от вида, структуры бетонного изделия и способа его тепловой обработки в теплотехнологических установках. В дополнении к граничным условиям при решении уравнения теплопроводности использовались модельные функции температуры от времени тепловой обработки, с помощью которых задавали температуру нагреваемой среды. В работе использовали табулированные функции температуры как модельные, так и экспериментально измеренные в различных областях бетонного изделия, при проведении экспериментальных исследований на предприятии «Завод КПД» ОАО «Строительно-монтажный трест № 16, г. Новополоцк».

В качестве основного уравнения для расчета распределения температуры в пространстве бетонного тела в течение времени ТВО при нагреве и гидратации бетона использовалось уравнение теплопроводности, учитывающее экзотермический характер реакции гидратации цемента, которое в указанной постановке имеет вид

дТ (х, у, г, т) дт

сьрь

Г д ^ Хь (Н (х, у, г, т), Т (х, у, г, т))дТ (^угТ) дх

д I 1ь (Н (х, у, г, т), Т (х, у, г, т))

дТ (х, у, г, т) ду

ду

(1)

д I 1ь (Н (х, у, г, т), Т (х, у, г, т))

дТ (х, у, г, т) дг

дг

+ РсР Ш(х, у, г,т), Т(х, у, г, т))

где т - время, с;

х, у, г - декартовы координаты (ось У направлена вертикально), м;

Т(х, у, г, т) - температура бетона в точке, имеющей координаты (х, у, г), в момент времени т, °С; рь - плотность бетона, кг/м3; рс - удельный расход цемента в бетоне, кг/м3; сь - удельная теплоемкость бетона, Дж/(кг-К);

Ш(х, у, г, т) - удельная теплота, выделившаяся при гидратации цемента в точке с координатами (х, у, г) на момент времени т;

Н(х, у, г, т) = Ш(х, у, г, т) / Штах - коэффициент (степень) гидратации цемента, где <2тах - удельная теплота полной гидратации цемента, Дж/кг;

1ь(Н, Т) - теплопроводность бетона, зависящая от степени его гидратации и температуры, Вт/(м-К); Р?(б, Т) - удельная мощность тепловыделения при гидратации цемента (Вт/кг), которая определяется по формуле

Р (Ш,Т) =

дШ(Т, тприв Ш, Т)) дт

(2)

где Ш - удельная теплота, выделившаяся при гидратации цемента к рассматриваемому моменту времени, Дж/кг;

1

+

+

+

+

Q(T, т) - функция тепловыделения цемента, Дж/кг;

Тприв(2, Т) - приведенное время гидратации цемента, с, которое определяется как функция обратная Q(T,т) при фиксированном значении Т(тприв(б, Т)= Q -1(Т, Q).

Более подробно основные положения, описывающие разработанную и используемую в данной статье теплофизическую модель, ее граничные условия, а также сведения о функциях тепловыделения и приведенном времени гидратации, изложены в ранее опубликованных работах автора с соавторами [19, 20].

На рисунке 1 на примере сечения модельного изделия кубической формы, подвергаемого ТВО на стенде горизонтального формования, показаны расчетные области модели, включающие пространство твердеющего бетона (^1), пространство опалубки (^2), пространство тепловой среды (^з), области нагревателей (й4). Каждая из перечисленных областей характеризовалась теплопроводностью, плотностью и удельной теплоемкостью материалов, включенных в расчет. Обрабатываемое изделие по бокам ограничено металлической формой, нижней гранью лежит на обогреваемом жидким теплоносителем столе, а сверху укрыто тепловлагоизолирующим съемным рулонным покрытием.

^2

Столешница

^2

Опалубка

^з

Паровоздушная среда

Нагреватели

Рисунок 1. - Схематичная область расчетного моделируемого пространства

Начальные условия для расчетной модели

ГТ(X, у, г, 0) = 20° С, для (х, у, г) ей1 и О2 и О3 х, у, г,0) = 0 Дж/кг , для (х, у, г) е Оц.

(3)

Конечно-разностные формулы для компьютерной модели. Сущность метода моделирования, использованного в работе, состоит в разбиении всего 3-мерного пространства модели на элементы прямоугольной (кубической) формы (рисунок 2). Каждый элемент имеет параметры, описывающие физическое состояние среды (теплоемкость, плотность, теплопроводность, температуру, степень гидратации). Все элементы расчетного пространства образуют взаимосвязанную пространственную ЗБ-систему, состояние и параметры которой пересчитываются на каждом шаге моделирования процесса нагрева и гидратации с учетом начальных и граничных условий. Каждый такт моделирования соответствует дискретному моменту времени, которое изменяется с шагом Дт. Состояние каждого элемента изменяется на каждом такте моделирования в зависимости от его текущего состояния и состояния соседних элементов на предыдущем такте моделирования. Формула, по которой вычисляется новое значение параметра элемента (в частности, температуры, коэффициента гидратации), является конечно-разностной аппроксимацией численного решения дифференциального уравнения теплопроводности.

На рисунке 2 выделены элементы расчетного пространства, использованные для расчета распределения температур и коэффициента гидратации на основе конечно-разностных аппроксимаций.

Так, для сеточных элементов группы (1, 2) температура и степень гидратации вычислялись по формулам (4), (5), для элементов группы (3) - по формуле (6), для элементов группы (4) - по формуле (7), для элементов группы (5) - по формуле (8), для элементов группы (6) - по формуле (9), для элементов группы (7) - по формуле (10).

Рисунок 2. - Схема разбиения моделируемого расчетного пространства на группы дискретных элементов Дх = Ду = Дг = 5 мм

На основе уравнения (1) была составлена конечно-разностная схема для численного расчета изменения температуры в дискретных элементах пространства твердеющего бетона. Для сеточных элементов группы (1, 2) (рисунок 2) использовалась следующая формула для расчета температуры:

Т(х, у, г, т + Дт) = Т(х, у, г, т) +

Ат

сьРь

Н(х + Ах, у, г, т) + Н(х, у, г, т) Т(х + Ах, у, г, т) + Т(х, у, г, т) ^ Т(х + Ах, у, г, т) - Т(х, у, г, т)

2

2

Ах

-1

Н(х-Ах, у,г,т) + Н(х, у, г,т) Т(х-Ах, у, г,т) + Т(х, у, г,т) | Т(х, у,г,т) -Т(х-Ах, у,г,т)

/Ах +

2 2 ) Ах

Н (х, у + Ау, г, т) + Н (х, у, г, т) Т(х, у + Ау, г, т) + Т(х, у, г, т) ^ Т(х, у + Ау, г, т) - Т(х, у, г, т)

-1,

Ч

2 2 ) Ау

Н(х, у - Ау, г, т) + Н(х, у, г, т) Т(х, у - Ау, г, т) + Т(х, у, г, т) ^ Т(х, у, г, т) - Т(х - Ау, у, г, т) „ 2 ' 2 ) Ау

Н(х, у, г + Аг, т) + Н(х, у, г, т) Т(х, у, г + Аг, т) + Т(х, у, г, т) ^ Т(х, у, г + Аг, т) - Т(х, у, г, т)

/Ау +

-1ь

^ 2 2 ) Аг

Н(х, у, г - Аг, т) + Н(х, у, г, т) Т(х, у, г - Аг, т) + Т(х, у, г, т) ^ Т(х, у, г, т) - Т(х, у, г - Аг, т)

2

2

Аг

/Аг + (4)

+ РсР Ш( х, у, г, т), Т (х, у, г, т))),

где Ат - шаг дискретизации по времени, с;

Ах, Ау, Аг - шаги дискретизации по пространству, м.

Остальные величины, входящие в уравнение (4), расшифрованы ранее.

Расчет распределения тепловыделений и соответствующих им коэффициентов (степени) гидратации цемента производился по формулам:

Ш( х, у, г, т + Ат) = Ш( х, у, г, т) + Рд (Ш( х, у, г, т), Т (х, у, г, т)) Ат Ш( х, у, г, т + Ат)

(5)

Н (х, у, г, т +Ат) = -

Шт

Расчет изменения температуры в дискретных элементах пространства стальной опалубки и конструкции оборудования, исключая граничные элементы (элемент 3), вычислялся по формуле

1

ь

1

+

+

T (х, у, z, t + Dt) = T (x, у, z, t) +

At-Kt

cst pst

T(x + Dx, у, z, t) - T(x, у, z, t) T(x, у, z, t) - T(x - Dx, у, z, t)

Dx

Dx

Dx+

/Ду +

/Dz 1,

(6)

T(x, у + Ду, z, t) - T(x, у, z, t) - T(x, у, z, t) - T(x - Ду, у, z, t)

Ду Ду .

T(x, у, z + Dz, t) - T(x, у, z, t) T(x, у, z, t) - T(x, у, z - Dz, t) Dz Dz _

где cst - удельная теплоемкость стали (475 Дж/кг-K); pst - плотность стали (7850 кг/м3); Ist - теплопроводность стали (44,5 Вт/(м-К)).

Расчет изменения температуры бетона на границе с греющей средой (Q - Q3) (верхняя горизонтальная поверхность бетона, элемент 4) вычислялся по формуле

T (x, у, z, t + Dt) = T (x, у, z, t) + Dt

cbpb

l ( H(x + Dx, у, z, t) + H(x, у, z, t) T(x + Dx, у, z, t) + T(x, у, z, t) 1 T(x + Dx, у, z, t) - T(x, у, z, t) _

b [ 2 ' 2 ) Dx

' H(x - Dx, у, z, t) + H(x, у, z, t) T(x - Dx, у, z, t) + T(x, у, z, t) ^ T(x, у, z, t) - T(x - Dx, у, z, t)

. 2 , 2 ) Dx

-1b

+ [aup (Tenv (t) - T(x У, Z, t)) -

' H(x, у - Ду, z, t) + H(x, у, z, t) T(x, у - Ду, z, t) + T(x, у, z, t) ^ T(x, у, z, t) - T(x - Ду, у, z, t)

Dx +

-1

•h , / , b 1 2 2 ) Ду

H(x, у, z + Dz, t) + H(x, у, z, t) T(x, у, z + Dz, t) + T(x, у, z, t) ^ T(x, у, z + Dz, t) - T(x, у, z, t)

/Ду+

2

2

Dz

-1 ( H(x, у, z - Dz, t) + H(x, у, z, t) T(x, у, z - Dz, t) + T(x, у, z, t) ^ T(x, у, z, t) - T(x, у, z - Dz, t) - b t 2 , 2 ) Dz

Dz + (7)

+ рсРс (0( х, у, г, х), Т (X, у, г, т))),

где (х, у, г), (х + Дх, у, г), (х - Дх, у, г), (х, у - Ду, г), (х, у, г + Дг), (х, у, г - Дг) - координаты середин конечно-разностных элементов в пространстве бетона;

(х, у + Ду, г) - координаты середины конечно-разностного элемента в пространстве паровоздушной среды;

Твпг(т) - функция изменения температуры паровоздушной среды от времени (Т(х, у + Ду, г) = Тет,(т)); аир - коэффициент теплоотдачи (амр = 5 Вт/(м2К)).

Расчет изменения температуры опалубки на границе с паровоздушной средой (Ф2 - ^3) с одной стороны и бетоном (Ф2 - ^3) - с другой для вертикальной боковой поверхности в направлении координатной оси X (элемент 5) проводился по формуле

T (x, у, z, t + Dt) = T (x, у, z, t) +

Dt

cst pst

asd (Tenv(t) -T(x,у,z,t))--

T(x, у, z, t) - T(x - Dx, у, z, t)

(1/ 1b (H(x-Dx, у, z,t),T(x-Dx, у, z,t)) +1/1st )Dx/2 T (x, у + Ду, z, t) - T (x, у, z, t) T (x, у, z, t) - T (x - Ду, у, z, t)

Dx +

Ду Ду

T(x, у, z + Dz, t) - T(x, у, z, t) T(x, у, z, t) - T(x, у, z - Dz, t)

Dz

Dz

/Ду+

/Dz 1,

+

X

+

b

+

где (х, у, г), (х - Ах, у, г), (х, у+Ау, г), (х, у - Ау, г), (х, у, г+Аг), (х, у, г-Аг) - координаты середин конечно-разностных элементов в пространстве опалубки;

(х + Ах, у, г) - координаты середины конечно-разностного элемента в пространстве паровоздушной среды;

Твпг(т) - функция изменения температуры паровоздушной среды от времени (Т(х+Ах, у, г) = Тет,(т)); а^ - коэффициент теплоотдачи (а^ = 7 Вт/(м2К)).

Расчет изменения температуры бетона на границе с опалубкой (^1 - ^2) для вертикальной боковой поверхности в направлении координатной оси X (элемент 6) проведен по формуле

Т (х, у, г, т + Ат) = Т (х, у, г, т) +

Ат

сьрь

Т (х + Ах, у, г, т) - Т (х, у, г, т)

(1/ 1ь (Н (х + Ах, у, г, т), Т (х + Ах, у, г, т)) +1/)Ах/2

-1ь

Н(х - Ах, у, г, т) + Н(х, у, г, т) Т(х - Ах, у, г, т) + Т(х, у, г, т) ^ Т(х, у, г, т) - Т(х - Ах, у, г, т)

2 2 ) Ах

Н (х, у + Ау, г, т) + Н (х, у, г, т) Т(х, у + Ау, г, т) + Т(х, у, г, т) ^ Т(х, у + Ау, г, т) - Т(х, у, г, *)

Ах +

-1ь

2 2 ) Ау

Н (х, у - Ау, г, т) + Н(х, у, г, т) Т(х, у - Ау, г, т) + Т(х, у, г, т) ^ Т(х, у, г, т) - Т(х - Ау, у, г, т) 2 '2 ) Ау

Н(х, у, г + Аг, т) + Н(х, у, г, т) Т(х, у, г + Аг, т) + Т(х, у, г, т) ^ Т(х, у, г + Аг, т) - Т(х, у, г, т)

/Ау +

-1ь

^ 2 2 ) Аг

Н(х, у, г - Аг, т) + Н(х, у, г, т) Т(х, у, г - Аг, т) + Т(х, у, г, т) ^ Т(х, у, г, т) - Т(х, у, г - Аг, т)

2

2

Аг

Аг +

+ РсРс Ш( х, у, г, т), Т (х, у, г, т))), (9)

где (х, у, г), (х - Ах, у, г), (х, у + Ау, г), (х, у - Ау, г), (х, у, г + Аг), (х, у, г - Аг) - координаты середин конечно-разностных элементов в пространстве бетона;

(х + Ах, у, г) - координаты середины конечно-разностного элемента в пространстве опалубки. Расчет изменения температуры стола (тот же материал, что и опалубка) на границе с нагревателем (^2 - ^4) (элемент 7) проведен по формуле

Т (х, у, г, т + Ат) = Т (х, у, г, т) +

Ат

Т(х + Ах, у, г, т) - Т(х, у, г, т) Т(х, у, г, т) - Т(х - Ах, у, г, т)

Ах

Ах

Ах +

- Т (х, у + Ау, г, т) - Т (х, у, г, т) , К* --—--—-- - (Т(х, у, г, т) - Т„(т))

Ау

/Ау +

(10)

+1

'Я*

Т(х, у, г + Аг, т) - Т(х, у, г, т) Т(х, у, г, т) - Т(х, у, г - Аг, т) Аг Аг

/Аг I,

где (х, у, г), (х + Ах, у, г), (х - Ах, у, г), (х, у + Ау, г), (х, у, г + Аг), (х, у, г - Аг) - координаты середин конечно-разностных элементов в пространстве опалубки;

(х, у - Ау, г) - координаты середины конечно-разностного элемента в пространстве нагревателя; ак - коэффициент теплоотдачи, соответствующий конвективному теплообмену между водяным теплоносителем и греющей стенкой (был принят ак = 250 Вт/(м2К)).

В конечно-разностной модели использовалось множество элементов, которые являются граничными одновременно по нескольким направлениям и с разными типами границ. Например, верхний угло-

+

1

+

+

1

+

вой элемент бетонного куба верхней гранью граничит с паровоздушной средой, а двумя боковыми гранями - с опалубкой. Конечно-разностная схема расчета этих элементов представляет собой соответствующие комбинации фрагментов итерационных уравнений, приведенных в [19, 20], и их аналогов по направлениям У и 2.

Описанные выше частные случаи конечно-разностных формул расчета температур для различных элементов модели были обобщены. Итоговая формула расчета изменения температуры в произвольных элементах пространства твердеющего бетона (^0, включая граничные элементы, имеет вид

Т (х, у, г, х + Ах) = Т (х, у, г, х) +

Л^ Г у+Ах-у-Ах +У+Ау "У-Ау + У+Аг-у-Аг +РсРс ^ у, г, т), Т(х, у, г, х)) 1, (11) сьрь У Ах Ау Аг )

где

У+Ах =

Г Н(х + Ах, у, г, х) + Н(х, у, г, х) Т(х + Ах, у, г, х) + Т(х, у, г, х) | Т(х + Ах, у, г, х) - Т(х, у, г, Г)

НI-

Ах

(х+Ах, у,г,х)еО1

Т (х + Ах, у, г, х) - Т (х, у, г, х)

У-Ах =

(1/Хъ ( Н (х, у, г, х),Т (х, у, г, х) ) +1/^ )Ах/2( х+Ах, у^

ГН(х-Ах,у,г,х) + Н(х,у,г,х) Т(х-Ах,у,г,х) + Т(х,у,г,х) |Т(х-Ах,у,г,х)-Т(х,у,г,х)

2

Т(х - Ах, у, г, х) - Т(х, у, г, х)

2

Ах

(х-Ах, у,г,х)еО1

У+Ау =

(1/ХЬ (Н (х, у, г, х),Т (х, у, г, х))+1/^ )Ах^х^, у^ . Г Н (х, у + Ау, г, х) + Н(х, у, г, х) Т(х, у + Ау, г, х) + Т(х, у, г, х) | Т(х, у + Ау, г, х) - Т(х, у, г, х)

2

2

Ау

(х, у+Ау,г,х)еО1

аир (Теш(Г) - Т(х, у, г, х))|

(х,у+Ау,г,х)еО3

У-Ау

1 Г Н (х, у -Ау, г, х) + Н (х, у, г, х) Т (х, у -Ау, г, х) + Т (х, у, г, х) 1Т (х, у, г, х) - Т (х, у -Ау, г, х)

ь У 2 ' 2 ) Ау

Т(х, у, г, х) - Т(х, у - Ау, г, х)

(х, у-Ау,г,х)еО1

У+Аг =

(1/ ХЬ (Н(х, у, г,х),Т(х, у, г,х)) +1/ Х„) Ау /2 (х,у-Ау,г,х)еа2 Г Н (х, у, г + Аг, х) + Н(х, у, г, х) Т(х, у, г + Аг, х) + Т(х, у, г, х) | Т(х, у, г + Аг, х) - Т(х, у, г, х)

^ У 2 ' '

Т (х, у, г + Аг, х) - Т (х, у, г, х)

2

Аг

(х, у,г+Аг,х)е О]

(1 / Хь (Н(х, у, г, х), Т(х, у, г, х)) +1/ Х^) Аг /2

У-Аг =

(х, у,г+Аг,х)еО2

Н(х, у, г - Аг, х) + Н(х, у, г, х) Т(х, у, г - Аг, х) + Т(х, у, г, х) | Т(х, у, г, х) - Т(х, у, г - Аг, Г)

2

Т(х, у, г, х) - Т(х, у, г - Аг, х)

(1/ Хъ (Н(х, у, г, х), Т(х, у, г, х)) +1/ ) Аг /2

2

(х, у,г-Дг,х)еО2

Аг

(х, у, г-Агд^О!

Итоговая формула расчета изменения температуры в элементах пространства опалубки (Ф2), включая граничные элементы, имеет вид

2

2

ъ

ъ

X

ъ

где

Т(х, у, г,т + Ат) = Т(х, у, г,т) + f ¡+^ + ¡+Ау ^ + ¡+Аг

ся*ря* V Ах

Ау

Аг

¡+Ах =

Т (х + Ах, у, г, т) - Т (х, у, г, т)

Ах

(х+Ах, у, гд)е^2

Т (х + Ах, у, г, т) - Т (х, у, г, т)

(1/ 1ь (Н(х + Ах, у, г, т), Т(х + Ах, у, г, т)) +1/ )Ах/2

(х+Ах, у, гд)еП1

(Тею (т) - Т (х, у, г, т) )|

(х+Ах,у^д^^ '

¡-Ах = 1

Т(х, у, г, т) - Т(х - Ах, у, г, т)

Ах (х-Ах, у^д^^

Т(х, у, г, т) - Т(х - Ах, у, г, т)

(1/ 1ь (Н (х -Ах, у, г, т),Т (х -Ах, у, г, т)) +1/ )Ах /2

С^ (Т ( х, у, г,т) - ТеШ (т))| (х-Ах, у^^

(х-Ах, у^д^^

¡+Ау = '

Т (х, у + Ау, г, т) - Т (х, у, г, т)

Ау (х, у+Ау, г,т)еИ2

Т (х, у + Ау, г, т) - Т (х, у, г, т)

(1/ 1ь (Н (х, у + Ау, г, т), Т (х, у + Ау, г, т)) +1/1 * )Ау /2 аир (Теп» (т) - Т ( х, у, г,т) )1 (х, у+Аy'г'T)eWз'

(х, у+Ау, гд)еИ1

¡-Ау = '

Т(х, у, г, т) - Т(х, у - Ау, г, т)

Ау

(х, у-Ау, гд)еИ2

Т(х, у, г, т) - Т(х, у - Ау, г, т)

(1/ 1ь (Н (х, у -Ау, г, т), Т (х, у -Ау, г, т)) +1/1)Ау/2

а№ (Т(х, у, г,т) - Т„(т))|(х,у_Ау,гд)бй4 ,

(х, у-Ау, г,т)еП1

¡+Аг = 1

Т (х, у, г + Аг, т) - Т (х, у, г, т)

Аг

(х, у, г+Агд^^

Т (х, у, г + Аг, т) - Т (х, у, г, т)

(1/ 1ь (Н(х, у, г + Аг,т),Т(х, у, г + Аг,т)) +1/ )Аг /2

(х, у,г+Агд)еП1

(Теп» (т) - Т (х, у, г, т) )|

(х,у,г+Агд)еП3 '

(12)

1

1

1

¡-Dz -

T (x, у, z, t) - T (x, у, z - Dz, t)

Dz

(x, y,z-Az,t)eW2

T(x, у, z, t) - T(x, у, z - Dz, t)

(1/ 1b (H(x, у, z-Dz,t),T(x, у, z-Dz,t)) +1/ 1st )Dz /2 asd (T(x, у, z,t) - Tenv(t))|(x,у,z-Az,t)eW3 '

(x, у, z-Dz,t)eWi

Алгоритм расчета процессов нагрева и гидратации. На основе вышеперечисленных формул была разработана САЕ (система автоматизированного моделирования) для расчета эволюции процесса гидратации бетонных изделий, написанная на языке C++, которая условно названа «ПроЭнергоБетон-1» (Процессы энергосбережения бетона).

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

В программе использован следующий укрупненный алгоритм моделирования процесса распределения температуры и степени гидратации:

Шаг 1. Разбить все моделируемое пространство (прямоугольное) на элементы с размерами Ax, Ду, Az, т.е. создать трехмерный массив Mt элементов с размерностью I (по Х), J (по Y), K (по Z).

Шаг 2. Задать начальные значения (формула (3)) параметров согласно соответствующим областям Q1, Q3, Q4 для момента времени т = 0.

Шаг 3. Создать новый массив Мт + Дт элементов с такой же размерностью, как Мт [I, J, К].

Шаг 4. Такт моделирования:

Шаг. 4.1. Проход по всем элементам массива MT[ij,k], где i = 0...I, j = 0...J, к = 0...К, и вычисление массива значений T, Q, Н по формулам (4)-(12) в зависимости от принадлежности элемента области Q1, Q2, Q3, Q4. Новые значения присваивать соответствующим элементам массива Мт+Дт [ij,k];

Шаг 4.2. Перенести полученные значения параметров из массива Мт + Дт [ij,k] в массив Мт[у,к] (или поменять массивы местами);

Шаг. 4.3. Итерационное изменение значения времени т = т +Дт;

Шаг 5. Повторять процесс моделирования (шаг 4) пока т не будет равно длительности цикла термообработки.

На рисунках 4-7 представлены фрагменты визуализации в трех проекциях численного моделирования распределений температуры и коэффициента (степени) гидратации в центральном сечении модельного кубического изделия, полученные с помощью разработанного алгоритма для различных интервалов времени, прошедших с момента начала ТВО.

Рисунок 4. - Фрагмент численного моделирования распределения температур в изделии размером 0,3x0,3x0,3 м в трех проекциях на момент времени тепловой обработки т = 4 ч

Рисунок 5. - Фрагмент численного моделирования распределения коэффициента гидратации в изделии размером 0,3x0,3x0,3м в трех проекциях на момент времени тепловой обработки т = 4 ч

Рисунок 6. - Фрагмент численного моделирования распределения температур в изделии размером 0,3x0,3x0,3м в трех проекциях на момент времени тепловой обработки т = 10 ч

Рисунок 7. - Фрагмент численного моделирования распределения коэффициента гидратации в изделии размером 0,3x0,3x0,3м в трех проекциях на момент времени тепловой обработки т = 10 ч

Заключение и выводы. На основании предложенной авторами теплофизической модели процесса гидратации бетонного изделия [19, 20], опирающейся на неравновесном уравнении теплопроводности, дополненном функцией источника гидратации, при различных вариантах начальных и граничных условий, характеризующих условия тепловой обработки в промышленных теплотехнологических установках, с учетом подходов, намеченных в [18], разработан алгоритм компьютерного моделирования процессов ускоренной гидратации.

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

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

ЛИТЕРАТУРА

1. Малинина, Л. А. Снижение энергетических затрат при производстве сборного железобетона за счет рационального выбора цементов, назначения эффективных режимов термообработки бетона и учета экзотермии / Л.А. Малинина // Тез. докл. Всесоюз. научно-практического семинара по экономии энергии при производстве сборных железобетонных конструкций и изделий / Госстрой СССР. - М., 1984. - С. 53-58.

2. Теоретико-практические аспекты эффективности добавок-ускорителей твердения бетона / Э.И. Ба-тяновский [и др.] // Перспективы развития новых технологий в строительстве и подготовке инженерных кадров : сб. науч. ст. / ГрГУ им. Я. Купалы. - Гродно, 2010. - С. 278-283.

3. Дворкин, Л.И. Критерий рационального использования тепловой энергии в производстве бетона и железобетонных изделий / Л.И. Дворкин, О.Л. Дворкин // Технология бетонов. - 2014. - № 2. -С. 32-35.

4. Бабицкий, В.В. Прогнозирование характеристик твердеющего тяжелого бетона / В.В. Бабицкий, С. Д. Семенюк, М.С. Бибик // Ресурсоекономнш матер1али, конструкцп, буд1вл1 та споруди : зб. наук. праць. - Р1вне, 2009. - Вип. 18. - С. 3-12.

5. Бибик, М.С. К возможности проектирования режима тепловой обработки бетона в ямных пропарочных камерах / М.С. Бибик, Н.В. Суходоева, В.В. Бабицкий // Строительная наука и техника. - 2009. -№ 2. - С. 58-63.

6. Бибик, М.С. Об энергосберегающих режимах тепловой обработки бетонных и железобетонных изделий / М.С. Бибик, В.В. Бабицкий // Строительная наука и техника. - 2010. - № 4. - С. 55-59.

7. Ушеров-Маршак, А.В. Информационная технология бетона ускоренного твердения / А.В. Ушеров-Маршак, А.Г. Синякин // Бетон и железобетон. - 1994. - № 6. - С. 2-4.

8. Лыков, А.В. Теория тепло- и массопереноса / А.В. Лыков, Ю.А. Михайлов. - М. ; Л. : Госэнергоиз-дат, 1963. - 536 с.

9. Дмитрович, А.Д. Тепло- и массообмен при твердении бетона в паровой среде / А.Д. Дмитрович. -М. : Стройиздат, 1967. - 243 с.

10. Марьямов, Н.Б. Тепловая обработка изделий на заводах сборного железобетона (процессы и установки) / Н.Б. Марьямов. - М. : Стройиздат, 1970. - 272 с.

11. Федосов, С.В. Применение методов математической физики для моделирования массо- и энергопереноса в технологических процессах строительной индустрии / С.В. Федосов, A.M. Ибрагимов, А.В. Гущин // Строительные материалы. - 2008. - № 4. - С. 65-67.

12. Ge, Zh. Predicting temperature and strength development of the field concrete : Rétrospective Theses and Dissertations [Electronic resource] / Zhi Ge // Iowa State University. - 2005. - Mode of access: https://lib.dr.iastate.edu/cgi/viewcontent. cgi?article=2729&context=rtd. - Date of access: 14.12.2018.

13. Моделирование прогрева стеновых панелей при термической обработке / С.В. Федосов [и др.] // Строительные материалы. - 2007. - № 2. - С. 86-87.

14. Аксенчик, К.В. Исследование тепло- и массообмена в бетонных плитах, подвергаемых тепловой обработке / К.В. Аксенчик, Н.И. Шестаков // Вестник Череповецкого государственного университета. - 2010. - № 4. - С. 63-67.

15. Z^koutsky, J. Effect of temperature on the early-stage hydration characteristics of Portland cement: A large-volume calorimetric study / J. Z^koutsky, V. TydlMt, R. Cherny // Construction and Building Materials. - 2012. - № 36. - P. 969-976.

16. Modelling of heat of hydration for thick concrete constructions - a note / Bennet Kuriakose [et al.] // Journal of Structural Engineering. - 2015. - Vol. 42, No. 4, October - November. - P. 348-357.

17. Фролов, С.В. Математическое моделирование процесса тепловлажностной обработки бетонных и железобетонных изделий / С.В. Фролов, А.В. Лагутин // Инженерно-физический журнал. - 2002. -Т. 75. - № 3.

18. Чичко, А.Н. Компьютерные системы моделирования физических процессов / А.Н. Чичко // Вестник БНТУ. - 2003. - № 2. - С. 42-48.

19. Повышение энергетической эффективности теплотехнологического оборудования для производства бетонных изделий на основе численного моделирования нестационарных процессов / А.М. Нияков-ский [и др.] // Энергетика. Известия высших учебных заведений и энерг. объединений СНГ. - 2019. -Т. 62. - С. 177-191.

20. Верификация нестационарной математической модели твердения бетона в теплотехнологических установках / А.М. Нияковский [и др.] // Наука и техника. - 2019. - Т. 18. - № 2. - С. 137-145.

Поступила 28.03.2019

ALGORITHM AND NUMERICAL SCHEME FOR THE SIMULATION OF UNSTEADY PROCESSES OF HEAT TREATMENT OF CONCRETE PRODUCTS

A. NIYAKOVSKI, Y. YATSKEVICH, A. CHYCHKO

On the basis of the developed thermophysical model a computer simulation algorithm of the concrete products' accelerated hydration process is proposed. The algorithm uses the heat equation in a three-dimensional formulation, taking into account the distributed source of internal heat emissions, and allows to perform calculation and visualization of non-stationary fields of temperature and the hydration coefficient during the hardening of concrete products in thermal technological installations. The algorithm is designed for the development of the concrete products' heat treatment energy-saving modes with complex 3D-geometry and heterogeneous composite structure.

Keywords: computer modeling, software package, mathematical modeling, kinetics of cement hydration, temperature fields, three-dimensional equation of thermal conductivity, non-stationary thermal conductivity, energy saving, heat technologies.