CC BY
30
УДК 536.24 doi.10.17238/issn 2227-6572.2016.2.126
ОЦЕНКА ЧИСЛЕННЫХ СВОЙСТВ СХЕМЫ ДИСКРЕТИЗАЦИИ ПРИ МОДЕЛИРОВАНИИ ТЕПЛОПРОВОДНОСТИ МЕТОДОМ ПЕРЕМЕННЫХ СОСТОЯНИЯ
В.А. Стенин*
*Северный (Арктический) федеральный университет имени М.В. Ломоносова
При математическом моделировании процессов теплообмена используют численные, точные и приближенные аналитические методы. Современные численные методы не вполне пригодны для применения в инженерной практике в силу сложности и неуниверсальности используемых алгоритмов и программ. Точные аналитические методы эффективны при решении только линейных дифференциальных уравнений. Приближенные методы, разрабатываемые для решения определенного класса задач, при исследовании ряда проблем неприемлемы. Поэтому совершенствование методов моделирования теплообмена, их комбинирование с привлечением математического аппарата смежных дисциплин и разработка новых методов, сопряженных с современными прикладными вычислительными комплексами и симуляторами, являются актуальными задачами. В технологии судостроения для моделирования используются математический аппарат теории автоматического управления и прикладные компьютерные программы (Matrix, Simulink, VisSim и др.). При составлении уравнений движения тепловых объектов применяемые методы в целом достоверно отражают физику процессов, однако не учитывают появляющиеся в переходных режимах нелинейности. Совместное использование метода переменных состояния (МПС) и вычислительных комплексов (в частности, VisSim) позволяет нивелировать нелинейность уравнений движения и улучшить качество процессов управления. В работе на основе МПС получено матричное уравнение, которое включает математическое описание процессов теплообмена и условия однозначности. Исследуется проблема оценки численных свойств схемы дискретизации при моделировании теплопроводности МПС. Предлагается методика преобразования обыкновенных дифференциальных уравнений к матричным уравнениям переменных состояния. С помощью уравнений переменных состояния и характеристических уравнений определяется устойчивость решений схемы дискретизации. На примере решения задачи нестационарной теплопроводности неограниченной пластины с граничными условиями первого рода проведена оценка аппроксимируемости и сходимости конечно-разностной схемы, представленной МПС.
Ключевые слова: разностная система, метод переменных состояния, устойчивость, сходимость, аппроксимация, точность решений.
Контактное лицо: Стенин Валерий Александрович, адрес: 164520, Архангельская обл., г. Северодвинск, ул. Капитана Воронина, д. 6; e-mail: v.stenin@narfu.ru
Методы математического моделирования тепловых процессов широко освещены в многочисленных публикациях [1, с. 74-231; 2, с. 102; 3, с. 31-194]. Однако при решении прикладных задач теплообмена возникают трудности практического применения стандартных методов в условиях действующего производства.
При математическом моделировании процессов теплообмена используют численные, точные и приближенные аналитические методы. Современные численные методы (конечных разностей, конечных элементов, граничных элементов и др. [1, с. 231; 2, с. 102; 3, с. 165]) не вполне пригодны для использования в инженерной практике в силу сложности и не универсальности используемых алгоритмов и программ. Точные аналитические методы (метод интегральных преобразований, метод функции Грина, тепловых потенциалов и др. [1, с. 74; 2, с. 102; 3, с. 31-116]) эффективны при решении только линейных дифференциальных уравнений. К настоящему времени известно значительное число разнообразных приближенных методов (методы Ритца, Треф-фтца, Канторовича, взвешенных невязок, кол-локаций, метод Галеркина и др. [2, с. 102; 3, с. 167-194]). Общим недостатком большинства из них является то, что они разрабатывались для решения определенного класса задач. К тому же для исследования массивов сложной формы, комплексных систем и многих других задач использование ряда методов неэффективно и даже неприемлемо. Поэтому совершенствование методов моделирования теплообмена, их комбинирование с привлечением математического аппарата смежных дисциплин и разработка новых методов, сопряженных с современными прикладными вычислительными комплексами и симулято-рами, являются актуальными задачами.
Характерным примером этого направления является работа И.В. и В.А. Кудиновых, где с помощью интегрального метода теплового баланса на основе введения фронта температурного возмущения при использовании до-
полнительных граничных условий получены аналитические решения задач теплопроводности с переменными начальными условиями, с переменными во времени граничными условиями и внутренними источниками теплоты [4, с. 5]. В статье С.В. Губарева и соавторов рассматривается математическая модель и численный метод решения дифференциальных уравнений в частных производных параболического типа на основе метода клеточных автоматов, как численного метода решения [5, с. 1]. В работе РФ. Марданова панельным методом рассчитывается уравнение Лапласа, когда результат отыскивается в виде суперпозиции фундаментальных решений [6, с. 5]. Для построения математического описания процесса Р.Ш. Мисбаховым и В.Е. Мизоновым применен метод ячеечного моделирования, использующий математический аппарат теории цепей Маркова [7, с. 1]. Основным оператором модели является матрица теплопроводности, преобразующая вектор теплового состояния через малые промежутки времени. А.Н. Ди-легинской и И.А. Данилушкиным предлагается методика решения задач теплопроводности операционным методом на основе структурной теории распределенных систем, позволяющая с помощью передаточных функций тепловых объектов с распределенными параметрами получить структуру сложного объекта, а также программная реализация полученного решения в среде программирования Matlab [8, с. 4-5]. В статье О.А. Тусниной приведен алгоритм численного решения дифференциального уравнения стационарной трехмерной теплопроводности [9, с. 91]. Дискретизация дифференциального уравнения проводилась методом контрольных объемов. Алгоритм решения дифференциального уравнения теплопроводности был реализован в вычислительном комплексе TEPL.
В данной работе предлагается моделировать процессы теплообмена в элементах судовых энергетических установок МПС совместно с программой-симулятором. В технологии судостроительного производства при проектиро-
вании и настройке систем регулирования энергетических установок широко используются для моделирования математический аппарат теории автоматического управления (в частности, частотный анализ, МПС) и прикладные компьютерные программы (Matrix, Simulink, VisSim и др.). При составлении уравнений движения тепловых объектов применяемые методы в целом достоверно отражают физику процессов, однако не учитывают появляющиеся в переходных режимах нелинейности. Совместное использование метода МПС и вычислительных комплексов (в частности, VisSim [10, c. 51]) позволяет нивелировать нелинейность уравнений движения и улучшить качество процессов управления. При исследовании МПС математических моделей реальных тепловых процессов непрерывные системы приближенно заменяются дискретными. Важной проблемой, возникающей при такой замене, является задача сохранения качественных характеристик исследуемых систем и численных свойств схем дискретизации (сходимости, апроксими-руемости, точности и устойчивости решений) [11, с. 137].
Рассмотрим численные свойства схем дискретизации на основе МПС. При использовании метода прямых нестационарное уравнение теплопроводности вида (в частности, одномерная задача) [12, с. 435]
ди д 2u — = a—7 dt дх2
(1)
где А(0 - основная матрица состояния системы; В(7) - матрица управления; С(7) - матрица связи переменных состояния с выходом системы; ^(0 - матрица входа; и(7) - вектор переменных состояния объекта; z(t) - вектор управляющих воздействий; у(7) - вектор выходных переменных.
Решение системы, удовлетворяющее начальным условиям и(0) = и для вектора состояния и выходного вектора имеет вид [14, с. 327]
t
и (t ) = ф(/}-м(0) + |ф(/-т) •Bz (x)-d т,
0
( t Л
y(t) = C Ф^)-и(0) + j^(t-т)•Bz(т)• dт +Dz(t),
где фундаментальная матрица Ф(7) = еА7 = = ехр А7.
Покажем порядок решения задачи (1) в системе матричных уравнений МПС. В методе прямых с конечно-разностной аппроксимацией для (1) разобьем интервал 0 < х < 5 узлами (0 < к < п) с шагом h. По формуле численного дифференцирования для внутренних линий 1 < / < п - 1 получим систему уравнений [12, с. 442]
dut dt
—'- = mut-l - 2mut + muj+l, 1 < t < n — 1. (3)
с начальными и краевыми условиями и(х, 0) = 0, 0 < х < 5; и(0, 7) = Ф(0, и(5 7) = Ф(0, приводится к системе дифференциальных уравнений первого порядка.
Для системы, описываемой совокупностью обыкновенных линейных дифференциальных уравнений, МПС может быть представлен следующими зависимостями [13, с. 343]:
J u(t) = ^(t) • u(t) + B(t) • z(t), I y(t) = C (t) • u(t) + D(t) • z(t),
(2)
Допустим, что дано распределение температуры по толщине неограниченной пластины в виде некоторой функции fx). Условия однозначности включают: теплопроводность материала пластины X = 45,4 Вт/(мК); плотность р = 7800 кг/м3; теплоемкость c = 502 Дж/кгК; толщину пластины 5 = 0,024 м; температуропроводность a = X/(pc) = 1,1595-10—5 м2/с. Начальное распределение температуры fx) = Тн = = 20 °С; температура ограничивающих поверхностей T0 = 600 °С.
В начальный момент времени ограничивающие поверхности мгновенно нагреваются до некоторой температуры T0, которая поддерживается постоянной на протяжении всего про-
Рис. 1. Расположение узлов в пластине (одномерная задача): 5 - толщина пластины; T- температура в узловой точке; h - расстояние между узлами
цесса нагрева. Пластина разбита, для примера, на п = 4 слоя (рис. 1).
Запишем систему дифференциальных уравнений в форме (3):
= mT0 - 2шT1 + mT2,
= mTl -2mT2 + ш^,
= mT2 - 2шT3 + mT0.
(4)
Представим систему уравнений (4) в виде переменных состояния (2):
+ [ш 0 m]•[T0 ],
T! -2ш ш 0
= ш -2ш ш
Tз 0 ш -2ш
'
Tз _
T2 =[0 1 0^
T
1Н
T
•[^0].
Здесь ш = X / (срк2), X, с, р, к - соответственно теплопроводность, теплоемкость, плотность, толщина слоя модели.
Для того, чтобы система (2) была устойчива, необходимо и достаточно, чтобы все соб-
ственные числа s. (/ = 1, ..., п) матрицы А имели отрицательные действительные части, т. е. лежали слева от мнимой оси плоскости комплексного переменного s. Собственные числа s. можно найти из характеристического уравнения [15, с. 293]
[А - sE] = 0, (5)
где Е - единичная матрица.
Точное решение задачи, с которым сравниваются результаты моделирования в МПС, находится методом разделения переменных в виде [16, с. 197]:
е( х ,Бо) =
T (X ,Бо) - ^
То - T
(6)
да
Е Ап Х C0S (Мп (1 - Х )) • еХР (-^2 •F0) ,
где Бо = ах/Я2 - число Фурье (безразмерное время); X = х/Я = 2х/5 - безразмерная толщина; Ап = (-1)п+1 ■ 2/ци, = (2п - 1) я/2 - коэффициенты ряда.
При численном моделировании рассматриваем одну половину пластины, т. к. задача является симметричной. Поскольку ряд, содержащийся в формуле (6), - сходящийся, то можно взять фиксированное количество членов ряда, обеспечивающих заданную точность. В расчетах, выполненных в Mathcad 14, принято число членов ряда, равное п = 1000. При решении задачи МПС пластина разбивалась на 10 и 100 слоев, что позволило оценить такие численные свойства схемы дискретизации, как сходимость и аппроксимируемость. Результаты численного решения нестационарной теплопроводности методом переменных состояния представлены на рис. 2, см. с. 130.
Оценка погрешности расчета нестационарной теплопроводности МПС в сравнении с точным решением приведена в таблице, см. с.130.
Максимальная относительная погрешность вычислений для Бо = 2 и п = 10 составляет примерно 4 %; с увеличением числа слоев (п = 100)
=1
0,0 0.2 0,4 0:6 0.8 1,0 1.2 1,4 1,6 1,8 2,0
Рис. 2. Численное решение нестационарной теплопроводности МПС: 0 - относительная температура; X - относительная координата
она снижается на порядок (0,4 %). В соответствии с Д. Ши [11, с. 154], дискретизованное уравнение аппроксимирует исходное дифференциальное уравнение, если погрешности стремятся к нулю при h ^ 0 (п ^ да), что под-
вой, и аппроксимирующей. Матричная оценка устойчивости показала, что система (2) по теореме Ляпунова устойчива.
Таким образом, проведенный анализ показал целесообразность использования МПС для реше-
погрешность расчета теплопроводности пластины методом переменных состояния, %
Fo п = 10 п = 100
X = 0,6 X = 1,0 X = 0,6 X = 1,0
0,1875 0,3703 0,6583 0,0046 0,0113
0,6875 0,5206 0,5537 0,0106 0,0037
1,0000 1,1312 1,1894 0,0343 0,0093
1,5000 2,3077 2,1807 0,1572 0,1590
2,0000 3,8961 3,1579 0,1350 0,4367
Примечание: п - число слоев модели.
тверждается таблицей. Согласно теореме Лак-са [11, с. 161], конечно-разностная схема обладает сходимостью в том и только в том случае, если схема является одновременно и устойчи-
ния задач теплообмена совместно с программным комплексом, так как МПС обеспечивает сохранение качественных характеристик исследуемых систем и численных свойств схем дискретизации.
Список литературы
1. Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М., 2003. 784 с.
2. Тимошпольский В.И., Трусова И.А. Менделев Д.В., Ратников П.Э. Анализ методов математического моделирования процессов теплообмена в промышленных печах для нагрева металла // Литье и металлургия.
2012. № 2. С. 102-107.
3. Швыдкий В.С., ЛадыгичевМ.Г., Шаврин В.С. Математические методы теплофизики. М., 2005. 232 с.
4. Кудинов И.В., Кудинов В.А. Аналитические решения параболических и гиперболических уравнений тепломассопереноса. М., 2013. 391 с.
5. Губарев С.В., Берг Д.Б., Добряк П.В. Математическая модель и численный метод для решения задач диффузии и теплопроводности // Соврем. проблемы науки и образования. 2013. № 6. С. 1-9.
6. Марданов Р.Ф. Численные методы решения плоской задачи теплопроводности. Казань, 2007. 23 с.
7. Мисбахов Р.Ш., Мизонов В.Е. Моделирование теплопроводности в составной области с фазовыми переходами // Вестн. Иванов. гос. энергет. ун-та. 2015. Вып. 4. С. 1-6.
8. Дилигенская А.Н., Данилушкин И.А. Математическое моделирование систем с распределенными параметрами. Самара, 2012. 65 с.
9. Туснина О.А. Теплотехнический расчет конструкций численными методами // Вестн. Моск. гос. строит. ун-та. 2013. № 11. С. 91-99.
10. ЛевинаГ.А., ЧесноковМ.А., Слепова С.В. Моделирование колебательных процессов в программе VisSim. Челябинск, 2010. 63 с.
11. Ши Д. Численные методы в задачах теплообмена. М., 1988. 544 с.
12. Боглаев Ю.П. Вычислительная математика и программирование. М., 1990. 544 с.
13. Деруссо П., Рой Р., Клоуз Ч. Пространство состояний в теории управления. М., 1970. 620 с.
14. Сигорский В.П. Математический аппарат инженера. Киев, 1975. 768 с.
15. Теория автоматического управления: в 2 ч. / под ред. А.А. Воронова. Ч. 1. Теория линейных систем автоматического управления. М., 1986. 367 с.
16. Пехович А.И., ЖидкихВ.М. Расчеты теплового режима твердых тел. Л., 1976. 351 с.
References
1. Samarskiy A.A., Vabishchevich P.N. Vychislitel'naya teploperedacha [Computational Heat Transfer]. Moscow, 2003. 784 p.
2. Timoshpol'skiy V.I., Trusova I.A. Mendelev D.V., Ratnikov P.E. Analiz metodov matematicheskogo modelirovaniya protsessov teploobmena v promyshlennykh pechakh dlya nagreva metalla [Analysis of Methods of Mathematical Modeling of Heat Transfer Processes in Industrial Furnaces for Metals Heating]. Lit'e i metallurgiya [Foundry and Metallurgy], 2012, no. 2, pp. 102-107.
3. Shvydkiy V.S., Ladygichev M.G., Shavrin V.S. Matematicheskie metody teplofiziki [Mathematical Methods of Thermal Physics]. Moscow, 2005. 232 p.
4. Kudinov I.V., Kudinov V.A. Analiticheskie resheniya parabolicheskikh i giperbolicheskikh uravneniy teplomassoperenosa [Analytical Solutions of Parabolic and Hyperbolic Equations of Heat and Mass Transfer]. Moscow,
2013. 391 p.
5. Gubarev S.V, Berg D.B., Dobryak P.V. Matematicheskaya model' i chislennyy metod dlya resheniya zadach diffuzii i teploprovodnosti [The Mathematical Model and Numerical Method for the Problems Solution of Diffusion and Heat Conduction]. Sovremennye problemy nauki i obrazovaniya [Modern Problems of Science and Education], 2013, no. 6, pp. 1-9.
6. Mardanov R.F. Chislennye metody resheniyaploskoy zadachi teploprovodnosti [Numerical Methods of Solution of the Plane Problem of Heat Conduction]. Kazan, 2007. 23 p.
7. Misbakhov R.Sh., Mizonov VE. Modelirovanie teploprovodnosti v sostavnoy oblasti s fazovymi perekhodami [Thermal Conductivity Modelling in the Composite Field with Phase Transitions]. Vestniklvanovskogo Gosudarstvennogo Energeticheskogo Universiteta [Vestnik of Ivanovo State Power Engineering University], 2015, no. 4, pp. 1-6.
8. Diligenskaya A.N., Danilushkin I.A. Matematicheskoe modelirovanie sistem s raspredelennymi parametrami [Mathematical Modeling of Systems with Distributed Parameters]. Samara, 2012. 65 p.
9. Tusnina O.A. Teplotekhnicheskiy raschet konstruktsiy chislennymi metodami [Thermotechnical Calculation of Structures by Numerical Methods]. VestnikMoskovskogo gosudarstvennogo stroitel'nogo universiteta [Proceedings of Moscow State University of Civil Engineering], 2013, no. 11, pp. 91-99.
10. Levina G.A., Chesnokov M.A., Slepova S.V. Modelirovanie kolebatel'nykh protsessov v programme VisSim [Simulation of Oscillatory Processes in the VisSim Program]. Chelyabinsk, 2010. 63 p.
11. Shi D. Chislennye metody vzadachakh teploobmena [Numerical Methods in Heat Transfer Problems]. Moscow, 1988. 544 p.
12. Boglaev Yu.P. Vychislitel'nayamatematikaiprogrammirovanie [Computational Mathematics and Programming]. Moscow, 1990. 544 p.
13. Derusso P., Roy R., Klouz Ch. Prostranstvo sostoyaniy v teorii upravleniya [State Space in the Control Theory]. Moscow, 1970. 620 p.
14. Sigorskiy VP. Matematicheskiy apparatinzhenera [The Mathematical Apparatus of an Engineer]. Kiev, 1975. 768 p.
15. Teoriya avtomaticheskogo upravleniya: v 2 ch. Ch. 1. Teoriya lineynykh sistem avtomaticheskogo upravleniya [Automatic Control Theory. Part 1. The Theory of Linear Systems of Automatic Control]. Ed. by A.A. Voronov. Moscow, 1986. 367 p.
16. Pekhovich A.I., Zhidkikh V.M. Raschety teplovogo rezhima tverdykh tel [Calculations of the Thermal Regime of Solids]. Leningrad, 1976. 351 p.
doi.10.17238/issn 2227-6572.2016.2.126
V.A. Stenin*
*Northern (Arctic) Federal University named after M.V. Lomonosov (Severodvinsk, Arkhangelsk region, Russian Federation)
EVALUATION OF NUMERICAL PROPERTIES OF THE SAMPLE CIRCUIT
IN THE HEAT CONDUCTION MODELING BY THE STATE-VARIABLE APPROACH
The numerical, accurate and approximate analytical methods are used in the mathematical modeling of heat transfer processes. Modern numerical methods are impractical in the engineering practice due to their complexity and nonuniversality of the algorithms and programs. Accurate analytical methods are effective only in the solution of linear differential equations. Approximate methods developed to solve a certain class of tasks in the study of some problems are unacceptable. Therefore, improving the heat transfer modeling techniques, their combining with the mathematical apparatus of interfacing disciplines and the development of new methods linked with the modern application computer systems and simulators are the urgent tasks. The mathematical apparatus of the automatic control theory and applied computer programs (Matrix, Simulink, VisSim, etc.) are used in the shipbuilding technology for modeling. When generation of equations of the motion of thermal objects the applied methods accurately reflect the physics of processes, but do not take into account the emerging non-linearity in the transition modes. Sharing of the state-variable approach (SVA) and computing systems (in particular, VisSim) allows us to level the non-linearity of the equations of motion and improve the quality of the control flow. In the work on the SVA basis the matrix equation is obtained; it includes the mathematical description of heat transfer processes and single-valuedness conditions. The problem of evaluation of numerical properties of the sample circuit in the heat conduction modeling of SVA is studied. The method of transformation of ordinary differential equations to the matrix equations of state variables is proposed. The solution stability of the sample circuit is defined by the equations of state variables and characteristic equations. In the context of solving a task of transient heat conduction of non-limiting plate with the boundary conditions of the first genus we evaluated the approximability and mesh scheme convergency provided by the state-variable approach.
Keywords: differential system, state-variable approach, stability, convergency, approximation, solution accuracy.
Received on January 29, 2016
Corresponding author: Valeriy Stenin, address: Captain Voronin str., 6, Severodvinsk, Arkhangelsk region, 164520, Russian Federation; e-mail: v.stenin@narfu.ru
