Численный алгоритм расчета температурных полей пневматических шин в процессе вулканизации Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

Фундаментальная и прикладная химия, химическая технология

УДК 611.53

Профессор С.Г. Тихомиров, доцент Ю.В. Пятаков,

(Воронеж. гос.ун-т инж. технол.) кафедра информационных и управляющих систем. тел. (4732) 55-38-75

профессор О.В. Карманова, доцент В.И. Молчанов

(Воронеж. гос.ун-т инж. технол.) кафедра химии и химической технолологии и органических соединений и переработки полимеров. тел. (4732) 49-92-37

Professor S.G. Tikhomirov, associate Professor Yu.V. Pyatakov

(Voronezh state university of engineering technology)

Department of information and control systems. phone (4732) 55-38-75

professor O.V. Karmanova, associate Professor V.I. Molchanov

(Voronezh state university of engineering technology) Department of chemistry and chemical technology, and organic compounds and polymer processing. phone (4732) 49-92-37

Численный алгоритм расчета температурных полей пневматических шин в процессе вулканизации

The numerical algorithm for calculating temperature fields of the pneumatic tires during vulcanization

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

Summary. In the article discussed the mathematical formulation and numerical algorithm for solving the problem of calculating the temperature field in the process vulcanizing of the product, whose the thermal characteristics are depended on the temperature. As a mathematical model considered the system of differential equations of heat conduction, taking into account the change in the coefficients of thermal conductivity and heat density in multilayer product of the temperature. The system of equations is solved for a given initial distribution of temperature and for a given (time-dependent) temperatures on the border of the product to the press-mold and to the diaphragm. On the border of the contacts of adjacent layers are given the condition of continuity of temperature and heat flux. Change of the thermal conductivity from the time is approximated by linear functions. The activation energy of the vulcanization process is determined on the basis of experimental data obtained in the control test samples using a reometer. Considering the function representing the corresponding integrals of the thermal conductivity, the original system of differential equations is transformed to an equivalent system of differential equations convenient for constructing numerical algorithms for solving the problem. The resulting system of partial differential equations derived using the method of finite-difference approximation is replaced by a system of algebraic equations. Solution of the system of algebraic equations is carried out under the scheme explicit difference approximation. In the article calculated the temperature field for the tire at given initial and boundary conditions. Stability and accuracy of the numerical algorithm for solving the problem is demonstrated by the calculations performed with different sampling step along the time and space coordinates. Assessment of the degree of completion of the process is carried out by calculated equivalent time for temperature values, accepted as equivalent. The developed algorithm is an important part of the algorithm for solving the problem of determining the optimal mode of curing, providing the required quality products at the lowest cost, which is very important in conditions of continuous growth in energy prices.

Ключевье слова: вулканизация, температурное поле, начальные и граничные условия, пневматическая шина, устойчивость решения.

Keywords: vulcanization, temperature field, initial and boundary conditions, the pneumatic tire, stability of the solution.

© Тихомиров С.Г., Пятаков Ю.В., Карманова О.В., Молчанов В.И., 2015

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

Произведем постановку задачи. В качестве математической модели процесса вулканизации рассмотрим систему уравнений теплового баланса вида:

скРкдТ(t,x)/dt = V (Т)VT(t, x)]+ qk(t,x), x e Vk, k=1,2,...,K; (1)

где ск, pk ,Xk (t) соответственно удельная теплоемкость, плотность, теплопроводность k-го слоя; T (t, x), qk (t, x) - температура и плотность тепловыделения в точке x k -го слоя в момент времени t: qk(t,x)= q?M • expE -(T(t,x)-T3)/(R • T(t,x)-T3)], qCyM - суммарное количество тепла, выделяемое в в k-м слое; Ek - энергия активации процесса вулканизации k-го слоя; R - универсальная газовая постоянная; Тэ - эквивалентная температура, к которой приводятся результаты неизотермической вулканизации [1-5].

Систему уравнений (1) дополним:

- начальными условиями:

Т (0, x) = То, (2)

где Т0 - начальная температура, которая зависит от времени года и соответствует температуре в сборочном цехе;

- граничными условиями:

Т(t,x)=p1 (t,x), x e S1; (3)

T(t, x) = p2 (t, x), x e S2, (4)

где S1 - поверхность контакта пресс-форма-теплоноситель, S2 - поверхность контакта диафрагма-теплоноситель;

- контактными условиями в точках, принадлежащим поверхностям соприкосновения разнородных слоев:

lim Т (t, x ') = lim Т (t, x"), (5)

x' ^x x" ^x

lim lk (T (t, x'))dT (t, x')/dn =

x'^x

= lim Д.+1 (T(t, x"))dT (t, x")/dn . (6)

В (5), (6) х' е ¥к , х"е¥к+1, х е Як+2, Sk+2 -поверхность контакта смежных слоёв Ук и Ук+\, п -вектор нормали к Sk+2; к=1,2,..., К-1; К - количество разнородных материалов, входящих в расчетную модель.

Зависимость теплопроводности от температуры будем аппроксимировать линейными функциями вида:

Л(т) = ак -Ьк•Т , (7)

где ак, Ьк - константы.

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

(Т2 + 273,15) • (Т1 + 273,15), t" -- •

Ek =-

Т - Т ±2 11

f

где I' - время достижения максимального значения динамического модуля М(0 при температуре Т\; I"- время достижения максимального значения динамического модуля М(0 при температуре Т2.

Степень завершенности процесса вулканизации X ) в точке х к -го слоя в момент времени t будем оценивать по формуле:

(8)

X tt ) = м (tэ (t))-M0

W Mmax - Mo '

где t3 (t) - эквивалентное время вулканизации при температуре 71=190 °С:

• (t )=[

El

R У 7\ +273,15 Т (г)+273,15 J^

- соответственно,

М(t) - текущее, М0, Мт минимальное и максимальное значение динамического модуля М при Т1=190 °С (рисуник 1), Т (г) - текущее значение температуры в точке

х в момент времени г в °С.

1200 1500 180(1

и с

Рисунок 1. Кинетика вулканизации к-го слоя по экспериментальным данным, полученным с помощью реометра. 1-график зависимости динамического модуля М(0 от времени при температуре Т1=190 °С; 2-график М(0 при температуре Т2=170 °С

e

Далее приведено решение системы уравнений (1)-(6). Введем в рассмотрение функции Л к (г), определенные следующим

1

образом: Лк (Т) = ^ Лк (г)йт .

Тогда система уравнений теплопроводности (1) примет вид:

скРкдТ(t, л)/ дt = А Лк (г, х) + qk (г, х), х еКк, (9)

где А = V2 - оператор Лапласа, Лк (г, х) = Лк (Т (t, х)). Начальное условие (2) примет вид:

Лк (0, х) = 0, к=1,2,...,К. (10)

Граничные условия (3), (4) примут вид: Л(г, х) = Л(1)(г, х), х е (11)

Л(/,х) = Л(2)(/,х), хе¿2, (12)

где Л(1)(г, х)=Лт (^(г, х)), т - номер слоя, граница

которого в точке х совпадает с ; Л(2>(г ,х)=Лп(<р2(г,х)), п - номер слоя, граница которого в точке х совпадает с Б2.

Контактные условия (5), (6) примут вид:

Лк1 (', х

х е Я,

= Л"'

х е

дЛ к (t, х)

дп

_дЛмМ

х е

дп

х е

(13)

(14)

Учитывая форму объекта моделирования, систему уравнений (9) удобно записать в цилиндрической системе координат (рисунок 2).

Рисунок 2. Объект моделирования

Пренебрегая в первом приближении формой рисунка протектора, можно считать, что шина имеет осесимметрическую форму. Тогда система уравнений (9) будет иметь вид:

дТ (г, г, к)_д2Лк ((, г, к) 1 + J

дí дг2 г дг дк2

скРг

' + qk(í, г, к)

х(г.

(г,

Решение системы уравнений (10)-( 15) будем осуществлять с помощью метода их конечно-разностной аппроксимации на системе точек |И, хи), где г(т) = г (т-1) + Аг; г (0) = 0; т = 1,2,...№,, Аг - шаг дискретизации по времени, х.. = (.,к..), 1 = 1,2,...Д; I = 1,2,...К. Точки хи, х2.,., х^ , ; /=1,2,.,У расположены на отрезках [ хи,хК^], хи е51, х , е52 (рисунок 3).

■V,

Г X X V - Г —■—1-1—г '

Н'Н» ±ь

!' 1 :

X

о

Рисунок 3. Расположение точек х]1 = (г|1, к..)

Выполним конечно-разностную аппроксимацию дифференциального оператора в правой части (15) в точке |(т), х..).

Для этого перенесем начало системы координат Огк в точку х.. и выполним ее ортогональное вращение таким образом, чтобы ось О'к' новой системы координат О'г'к' совпала с прямой, проходящей через точки х.+и, х..,

х!-и (рисунок 4).

Введем обозначения:

Лг, = дЛ / дг ', Лг = д2Л / дг '2, Лк, = дЛ / дк', ЛИИ =д2Л / дк'2, Л= д2Л/дг'/дк'.

Рисунок 4. Расположение точек в системе координат О 'г 'к'.

где к = х1.

0

Разложим функцию Лк (г, х) в окрестности точки xj , в ряд Тейлора и, ограничившись

членами второго порядка малости относительно Ак", Ак2, получим систему уравнений для определения Лк, Лкк :

|Ак1 •Лк+Ак12/2 •Лкк= /l, (16) [Ак2 -Лк+Ак22/2 кк = /2, где Ак1 = к-1,,- - к", , Ак2 = к"+1,, - к",,, /1 =Лк (г(т), х-1, )-Лк (г(т), х,),

/2 =Лк (г(т), х+1, )-Лк (г(т), х,).

Аналогично, используя соответствующие разложения Лк (г(т), х), получим систему уравнений для определения Лг,, ЛгУ, Лгк:

Ат3 •ЛГ,+Агз •Акз • Л^ +Ат3 /2• Л„ = /, АГ4 • Л г +АГ4 • Ак А • Л гу +Аг„2 /2 • Л „ = /4; Аг5 •Лг, +АТ5 • Ак5 • Л,, +Аг52 /2• Л„ = /,

(17)

где Ак3 = к-Ц+1 -к",, , Аг3 = Г-1,,+1 - О", , Ак4 = к'и+1 -к",, Аг4 = г;,,+1- г",,, Ак5 = к,,-- к",, Аг5=г;,,-1- ;,

/з =Лк ((т), х;-1,+1 )-Лк (г(т), х„ )-АкзЛк,-Акз2/2-Акк,

/4 =Лк (/"', х;,+1 )-Лк (/"', )-Ак4 •Лк-Ак42/2 .Лкк,

/5 =Лк(г(т),х;,-1 )-Лк(г(т),)-Ак5 •Лк -Ак52/2•Лкк . Учитывая инвариантность оператора

д2Лк (í, г, к) д2Лк (г, г, к)

——2—'-+——г—'- относительно переноса

дг дк

начала системы координат Огк и относительно ее ортогонального вращения, вокруг точки О' получим:

д2Лк (г, г, к) д2Л к (г, г, к)_д2Лк (г, Г, к ') д2Лк (г, г', к ')

дг2 дк2 дг'2

Соответственно: дЛк (г, г, к)/ дг = Лг • п + Лк, • п2,

дк

-. (18)

(19)

где п1, п2 - направляющие косинусы вектора

Ог в системе координат О 'г 'к'.

Таким образом, выражение для аппроксимации дифференциального оператора в правой части (15) будет иметь вид:

ахкк)+1 Мтк)+9к ьА ( «х..)+

я,.2 я,. яг,2 «\> ' > ;,1 к\ щ

+Ь;,+1Лк (г(т), х;,+1)+ь.-Л (г«, х)+Ь;_1^+1Лк (г(т), х ^)+ + Ь;-и Лк (г ^ х;-1,, )+ Ь;+1,, Лк (гЧ х;+1, )+ Чк (гЧ х;,, ) где коэффициенты Ь;, определяются подстановкой вычисленных, путем решения систем уравнений (16), (17), значений Лк,, Лкк , Лг,, ЛгУ, Лгк в выражения (18), (19).

Учитывая, что Хк (г) является медленно меняющейся функцией температуры, решение системы (10)-( 15) во внутренних точках к-го слоя (х;, е¥к \дКк ) при шаге дискретизации по времени:

А < екРк А2/(4ак ) (20)

можно осуществить, используя явную разностную схему:

Т(г(т+1),х)= Т(гм,хи )+Аг(скРкЛк (г(т),хи)+ + Ь;, + Л (г(т), х;,+ 1 )+ Ь;,_Л х)+ Ь^+Л (г^ х^ ) + +Ь;_иЛк (г(т), х;_и)+Ь;+иЛк (г(т), хт)+ Чк (г(т), хи)], (21) где А - минимальное расстояние между точками х;,; ак -параметры аппроксимации теплопроводности к-го слоя. Для точек, принадлежащих поверхностям смежных слоев (х;, е Sk+2),

аппроксимацию дифференциального оператора в (15) будем осуществлять по аналогичной схеме на основании метода баланса [6].

Приведем пример расчета. В примере выполнен расчет температурного поля для шины (рисунок 5), имеющей наружный диаметр -1,078 м и посадочный диаметр - 0,6096 м.

Значения удельной плотности и теплоемкости:

- для протектора р=1127; с1=717;

- для брекерар2=1106; с2=641;

- для каркаса р3=1071; с3=615;

- для боковины р4=1134; с4=663;

- для борта р5 =7700; с5 =460;

- для наполнительного шнура р6 =1118; с6=664.

Размерность величин рк, к=1,2,_,5 дана в

кг/м3. Размерность величин ск, к=1,2,_,5 дана

в Дж/(кг град).

Рисунок 5. Конструкция пневматической шины. 1 - протектор; 2 - брекер; 3 - каркас; 4 - боковина; 5 - борт; 6 - наполнительный шнур

Параметры аппроксимации коэффициентов теплопроводности (7):

- для протектора а1 =0,1612; Ь1 =0,0002;

- для брекера а2=0,1793; Ь2=0,0003;

- для каркаса а3 =0,07; Ь3 =0,0;

- для боковины а4 =0,1941; Ь4 =0,0004;

- для борта а5=50,0; Ь5=0,0;

- для наполнительного шнура а6 =0,3426; Ь6 =0,0006.

Значения коэффициентов ак, к=1,2,...,5 даны в Вт/(мград). Значения коэффициентов Ьк, к=1,2,...,5 даны в Вт/(мградс).

При заданных начальных: Т0 =20 °С и граничных:

[150 °С, при 0 < t < 2700 с, 20 °С, при t > 2700 с условиях получено решение в виде изменения температурного поля во времени.

Расчеты температуры проводились в точках хп , / = 1,2,...Д; I = 1,2,...К (рисунок 3) при

N=58; К=15, /=1,2,...,20; К21=13; К22=11; К23=9; К24=7; К25=7; К26=6; К=5, 1=27, 28 , . ,49; К50=7; К51=8; К=9, 1=52, 53 , ... ,57; К58=1. Общее количество точек хл равно 557.

В таблице 1 приведены результаты расчетов, выполненных с разными значениями шага дискретизации по времени Аг. Расчеты проводились в точках х. 1, 1=1,2,.. .,15.

Т а б л и ц а 1

Результаты расчетов температурного поля

в момент времени г=1800 с, выполненные с разным значением Аг

<р1 (г, х) = р2 (г, х) =

Номер точки Значения шага дискретизации по времени Аг , с

8 4 2 1

1 150,0 150,0 150,0 150,0

2 149,4 149,5 149,5 149,5

3 149,0 149,0 149 ,0 149,0

4 148,6 148,5 148,5 148,5

5 148,2 148,2 148,1 148,1

6 147,9 147,9 147,9 147,9

7 147,8 147,7 147,7 147,7

8 147,7 147,7 147,6 147,6

9 147,8 147,7 147,7 147,7

10 147,8 147,8 147,8 147,8

11 148,0 148,0 148,0 148,0

12 148 ,4 148,4 148,4 148,4

13 148,9 148,9 148,9 148,9

14 149,4 149 ,4 149,4 149,4

15 150,0 150,0 150,0 150,0

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

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

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

отрезках [ х11., хК^ ] последовательно увеличивалось в 2, 4 и 8 раз. Величина А при этом уменьшалась соответственно 2, 4 и 8 раз.

Шаг дискретизации Аг определялся в соответствии с условием (20). Значения температуры в таблице 2 рассчитаны для тех же точек, что и в таблице 1.

Т а б л и ц а 2 Результаты расчетов температурного поля в момент времени г= 1800 с, выполненные с разным значением А

Номер точки Значения А, м

1,6х10"3 8,0х10"4 4,0х10"4 2,0х10"4

1 150,0 150,0 150,0 150,0

2 149,5 149,2 149,0 148,9

3 149,0 148,4 148,1 148,0

4 148,5 147,7 147,3 147,1

5 148,1 147,7 146,6 146,3

6 147,9 146,7 146,0 145,7

7 147,7 146,4 145,7 145,2

8 147,6 146,3 145,5 145,0

9 147,7 146,3 145,4 145,0

10 147,8 146,4 145,6 145,2

11 148,0 146,8 146,1 145,7

12 148,4 147,3 146,7 146,4

13 148,9 148,1 147,6 147,4

14 149,4 148,9 148,7 148,6

15 150,0 150,0 150,0 150,0

Примечание. В таблице 2 подчеркнуты расхождения с результатами, выполненными с шагом дискретизации А =2,0х10"4 м.

Приведенные в таблице2 результаты показывают устойчивую сходимость решения с уменьшением величины А .

На рисунках 6-10 приведены распределения температуры и степени завершенности процесса X(г) (см. формулу (8)) в шине для различных значений времени.

Рисунок 6. Температура (в °С) в момент времени г=1800 с

I I I I

14 \\\ \\\

\\

Рисунок 7. Степень завершенности процесса X (г) в момент времени /=1800 с

Рисунок 8. Температура (в °С) в момент времени г=3000 с

На момент времени г=3600 с степень завершенности процесса вулканизации во всех точках шины практически равна 1.

Заключение:

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

2) рассмотрен способ конечно-разностной аппроксимации системы диффе-

в момент времени 1=3000 с

г=3600 с

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

3) выполнены расчеты температурного поля для шины при заданных начальных и граничных значениях температуры;

4) получены результаты, демонстрирующие устойчивость и точность решения;

5) для различных моментов времени г выполнен расчет степени завершенности процесса вулканизации в изделии.

ВестникВГУИТ, №2, 205

ЛИТЕРАТУРА

REFERENCES

1 Лукомская А.И., Минаев Н.Т., Кепер-ша Л.М., Милкова Е.М. Оценка степени вулканизации резин в изделиях. М.: ЦНИИТ-Энефтехим, 1972. 43 с.

2 Денисов А. П., Громов Ю. Ю., Покорный Ю. В. Исследование процесса вулканизации при местном ремонте шин. // ВГУ, МГУ. Воронеж. Информационные технологии и системы. 1999. С. 97.

3 Карманова О. В., Тихомиров С. Г., Пятаков Ю.В., Касперович А. В. и др. Моделирование кинетики неизометрической вулканизации массивных резиновых изделий // Труды БГТУ. 2014. № 4. С. 100-104.

4 Ищенко В. А., Шаптала М. В. Особенности расчетов режимов вулканизации пневматических шин с учетом трехмерности конструкции // Системные технологии: региональный межвуз. сб. науч. трудов. 2008. Вып. 2 (55). С. 147 - 158.

5 Власко А.В., Сахаров М.Э., Пориц-кая З. Влияние неизотермической вулканизации и механические свойства резиновых и ре-зинокордных образцов // Каучук и резина. 1991. № 6. С. 6-8.

6 Самарский А. А., Тихонов А. Н. Уравнения математической физики : учебное пособие. М.: Изд-во МГУ, 1999. 798 с.

1 Lukomskaya A.I., Minaev N.T., Kepersha L.M., Milkova E. M. Otsenka stepeni vulkanizatsii rezin v izdeliyakh [Assessment of the degree of vulcanization of rubber in products]. Moscow, TsNII I Enefekhim, 1972. 43 p. (In Russ.).

2 Denisov A. P., Gromov Yu. Yu., Pokornyi Yu. V. Investigation of the vulcanization process at the local tire repair VGU, MGU. Voronezh. Informatsionnye tekhnologii i sistemy. [Voronezh State University, Moscow State University. Voronezh. Information technologies and systems], 1999, pp. 97. (In Russ.).

3 Karmanova O. V., Tikhomirov S. G., Py-atakov Yu.V., Kasperovich A. V. et al. Modelling of the kinetics of non-isometric massive rubber vulcanization. Proceedings BSTU. [Trudy BGTU], 2014, no. 4, pp. 100-104. (In Russ.).

4 Ishchenko V. A., Shaptala M. V. Features of the operating conditions of vulcanization of pneumatic tires with the three-dimensional structure. Sis-temnye tekhnologii: regional'nyi mezhvuz. sb. nauch. trudov. [System Technology: Regional Interuniversity collection of scientific papers], 2008, issue 2 (55), pp. 147 - 158. (In Russ.).

5 Vlasko A. V., Sakharov M. E., Poritskaya Z. Effect of non-isothermal curing and mechanical properties of rubber and rubber-cord samples. Kauchuk i rezi-na. [Rubber and rubber], 1991, no. 6, pp. 6-8. (In Russ.).

6 Samarskii A.A., Tikhonov A.N. Uravneniya matematicheskoi fiziki [Equations of mathematical physics]. Moscow, Izd-vo MGU, 1999. 798 p. (In Russ.).