All issues

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

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

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

This paper presents algorithm for modeling set of trajectories of non-stationary time series, based on a numerical scheme for approximating the sample density of the distribution function in a problem with fixed ends, when the initial distribution for a given number of steps transforms into a certain final distribution, so that at each step the semigroup property of solving the Liouville equation is satisfied. The model makes it possible to numerically construct evolving densities of distribution functions during random switching of states of the system generating the original time series.

The main problem is related to the fact that with the numerical implementation of the left-hand differential derivative in time, the solution becomes unstable, but such approach corresponds to the modeling of evolution. An integrative approach is used while choosing implicit stable schemes with “going into the future”, this does not match the semigroup property at each step. If, on the other hand, some real process is being modeled, in which goal-setting presumably takes place, then it is desirable to use schemes that generate a model of the transition process. Such model is used in the future in order to build a predictor of the disorder, which will allow you to determine exactly what state the process under study is going into, before the process really went into it. The model described in the article can be used as a tool for modeling real non-stationary time series.

Steps of the modeling scheme are described further. Fragments corresponding to certain states are selected from a given time series, for example, trends with specified slope angles and variances. Reference distributions of states are compiled from these fragments. Then the empirical distributions of the duration of the system’s stay in the specified states and the duration of the transition time from state to state are determined. In accordance with these empirical distributions, a probabilistic model of the disorder is constructed and the corresponding trajectories of the time series are modeled.

Проведен сравнительный анализ двух численных методик моделирования нестационарных режимов термогравитационной конвекции и теплового поверхностного излучения в замкнутой дифференциально обогреваемой кубической полости. Рассматриваемая область решения имела две изотермические противоположные вертикальные грани, остальные стенки являлись адиабатическими. Поверхности стенок считались диффузно-серыми, т. е. их направленные спектральные степень черноты и поглощательная способность не зависят ни от угла, ни от длины волны, но могут зависеть от температуры поверхности. Относительно отраженного излучения использовались два предположения: 1) отраженное излучение является диффузным, т. е. интенсивность отраженного излучения в любой точке границы поверхности равномерно распределена по всем направлениям; 2) отраженное излучение равномерно распределено по каждой поверхности замкнутой области решения. Математическая модель, сформулированная как в естественных переменных «скорость–давление», так и в преобразованных переменных «векторный потенциал–вектор завихренности», реализована численно методом контрольного объема и методом конечных разностей соответственно. Следует отметить, что анализ радиационного теплообмена проведен с использованием метода сальдо в варианте Поляка.

При решении краевой задачи в естественных переменных методом контрольного объема для аппроксимации конвективных слагаемых применялся степенной закон, для диффузионных слагаемых — центральные разности. Разностные уравнения движения и энергии разрешались на основе итерационного метода переменных направлений. Для поиска поля давления, согласованного с полем скорости, применялась процедура SIMPLE.

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

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

Comparative analysis of two numerical methods for simulation of unsteady natural convection and thermal surface radiation within a differentially heated cubical cavity has been carried out. The considered domain of interest had two isothermal opposite vertical faces, while other walls are adiabatic. The walls surfaces were diffuse and gray, namely, their directional spectral emissivity and absorptance do not depend on direction or wavelength but can depend on surface temperature. For the reflected radiation we had two approaches such as: 1) the reflected radiation is diffuse, namely, an intensity of the reflected radiation in any point of the surface is uniform for all directions; 2) the reflected radiation is uniform for each surface of the considered enclosure. Mathematical models formulated both in primitive variables “velocity–pressure” and in transformed variables “vector potential functions – vorticity vector” have been performed numerically using finite volume method and finite difference methods, respectively. It should be noted that radiative heat transfer has been analyzed using the net-radiation method in Poljak approach.

Using primitive variables and finite volume method for the considered boundary-value problem we applied power-law for an approximation of convective terms and central differences for an approximation of diffusive terms. The difference motion and energy equations have been solved using iterative method of alternating directions. Definition of the pressure field associated with velocity field has been performed using SIMPLE procedure.

Using transformed variables and finite difference method for the considered boundary-value problem we applied monotonic Samarsky scheme for convective terms and central differences for diffusive terms. Parabolic equations have been solved using locally one-dimensional Samarsky scheme. Discretization of elliptic equations for vector potential functions has been conducted using symmetric approximation of the second-order derivatives. Obtained difference equation has been solved by successive over-relaxation method. Optimal value of the relaxation parameter has been found on the basis of computational experiments.

As a result we have found the similar distributions of velocity and temperature in the case of these two approaches for different values of Rayleigh number, that illustrates an operability of the used techniques. The efficiency of transformed variables with finite difference method for unsteady problems has been shown.

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

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

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

Сформулированная в работе математическая модель позволяет описывать процесс распространения пламени по поверхности горючего материала. Газодинамика моделируется системой уравнений Навье – Стокса, вихревое течение описывается комбинированной моделью турбулентности RANS–LES (DDES), турбулентное горение — комбинированной моделью горения Eddy Break-Up с учетом кинетических эффектов, теплопередача излучением — методом сферических гармоник первого порядка аппроксимации (P1). Решение уравнений производится в программном пакете OpenFOAM.

Reduction of the fire hazard of polymeric materials is one of the important scientific and technical problems. Since complexity of experimental procedures associated with flame spread, establishing reacting flows theoretical basics turned out to be crucial field of modern fundamental science. In order to determine parameters of flame spread over solid combustible materials numerical modelling methods have to be improved. Large amount of physical and chemical processes taking place needed to be resolved not just separately one by one but in connection with each other in gas and solid phases.

Upward flame spread over vertical solid combustible material is followed by unsteady eddy structures of gas flow in the vicinity of flame zone caused by thermal instability and natural convection forces accelerating hot combustion products. At every moment different amount of heat energy is transferred from hot gas-phase flame to solid material because of eddy flow structures. Therefore, satisfactory heat flux and eddy flow modelling are important to estimate flame spread rate.

In the current study we evaluated parameters of numerical method for flame spread over solid combustible material problem taking into account coupled nature of complex interaction between gas phase, solid material and eddy flow resulted from natural convection. We studied aspects of different approximation schemes used in differential equations integration process over space and time, of fields relaxation during iterations procedure carried out inside time step, of different time step values.

Mathematical model formulated allows to simulate flame spread over solid combustible material. Fluid dynamics is modeled by Navier – Stokes system of equations, eddy flow is described by combined turbulent model RANS–LES (DDES), turbulent combustion is resolved by modified turbulent combustion model Eddy Break-Up taking into account kinetic effects, radiation transfer is modeled by spherical harmonics method of first order approximation (P1). The equations presented are solved in OpenFOAM software.

В статье представлены результаты численного моделирования вихревого пространственного нестационарного движения среды, возникающего около боковой и донной поверхностей десантного модуля при его спуске в атмосфере Марса. Численное исследование проведено для высокоскоростного режима обтекания при различных углах атаки. Математическое моделирование осуществлено на основе модели Навье – Стокса и модели равновесных химических реакций для газового состава марсианской атмосферы. Результаты моделирования показали, что при рассматриваемых условиях движения спускаемого аппарата около его боковой и донной поверхностей реализуется нестационарное течение, имеющее ярко выраженный вихревой характер. Численные расчеты указывают на то, что в зависимости от угла атаки нестационарность и вихревой характер потока могут проявляться как на всей боковой и донной поверхностях аппарата, так и, частично, на их подветренной стороне. Для различных углов атаки приводятся картины вихревой структуры потока около поверхности спускаемого аппарата и в его ближнем следе, а также картины полей температуры и показателя адиабаты. Нестационарный характер обтекания подтверждается представленными временными зависимостями газодинамических параметров потока в различных точках поверхности аппарата. Проведенные параметрические расчеты позволили построить зависимости аэродинамических характеристик спускаемого аппарата от угла атаки. Математическое моделирование осуществляется на основе являющегося методом конечных объемов консервативного численного метода потоков, основанного на конечно-разностной записи законов сохранения аддитивных характеристик среды с использованием upwind-аппроксимаций потоковых переменных. Для моделирования возникающей при обтекании сложной вихревой структуры потока около спускаемого аппарата используются неравномерные вычислительные сетки, включающие до 30 миллионов конечных объемов с экспоненциальным сгущением к поверхности, что позволило выявить мелкомасштабные вихревые образования. Численные исследования проведены на базе разработанного комплекса программ, основанного на параллельных алгоритмах используемого численного метода и реализованного на современных многопроцессорных вычислительных системах. Приведенные в статье результаты численного моделирования получены при использовании до двух тысяч вычислительных ядер многопроцессорного комплекса.

The article presents the results of numerical modeling of the vortex spatial non-stationary motion of the medium arising near the lateral and bottom surfaces of the descent module during its movement in the atmosphere of Mars. The numerical study was performed for the high-speed streamline regime at various angles of attack. Mathematical modeling was carried out on the basis of the Navier – Stokes model and the model of equilibrium chemical reactions for the Martian atmosphere gas. The simulation results showed that under the considered conditions of the descent module motion, a non-stationary flow with a pronounced vortex character is realized near its lateral and bottom surfaces. Numerical calculations indicate that, depending on the angle of attack, the nonstationarity and vortex nature of the flow can manifest itself both on the entire lateral and bottom surfaces of the module, and, partially, on their leeward side. For various angles of attack, pictures of the vortex structure of the flow near the surface of the descent vehicle and in its near wake are presented, as well as pictures of the gas-dynamic parameters fields. The non-stationary nature of the flow is confirmed by the presented time dependences of the gas-dynamic parameters of the flow at various points on the module surface. The carried out parametric calculations made it possible to determine the dependence of the aerodynamic characteristics of the descent module on the angle of attack. Mathematical modeling is carried out on the basis of the conservative numerical method of fluxes, which is a finitevolume method based on a finite-difference writing of the conservation laws of additive characteristics of the medium using «upwind» approximations of stream variables. To simulate the complex vortex structure of the flow over descent module, the nonuniform computational grids are used, including up to 30 million finite volumes with exponential thickening to the surface, which made it possible to reveal small-scale vortex formations. Numerical investigations were carried out on the basis of the developed software package based on parallel algorithms of the used numerical method and implemented on modern multiprocessor computer systems. The results of numerical simulation presented in the article were obtained using up to two thousand computing cores of a multiprocessor complex.

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

Построена математическая модель обтекания затупленного тела сверхзвуковым потоком вязкого газа с твердыми частицами. Решение системы нестационарных уравнений Навье–Стокса в консервативных переменных осуществляется бессеточным методом, в основе которого лежит аппроксимация частных пространственных производных газодинамических величин и содержащих их функций методом наименьших квадратов на множестве распределенных в области расчета узлов. Расчет невязких потоков выполняется методом HLLC в сочетании с MUSCL-реконструкцией третьего порядка, вязких потоков — схемой второго порядка. МНК-аппроксимация частных производных параметров газа по направлению также применяется для реализации краевых условий Неймана на выходной границе области расчета, а также поверхностях обтекаемых тел, которые считаются изотермическими твердыми стенками.

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

Выполнена программная реализация представленных алгоритмов на основе технологии параллельных гетерогенных вычислений OpenCL. Представлены результаты моделирования движения крупной частицы вдоль оси симметрии сферы навстречу набегающему потоку с числом Маха $\mathrm{M}=6$.

The work is devoted to numerical modeling of two-phase flows, namely, the calculation of supersonic flow around a blunt body by a viscous gas flow with an admixture of large high inertia particles. The system of unsteady Navier – Stokes equations is numerically solved by the meshless method. It uses the cloud of points in space to represent the fields of gas parameters. The spatial derivatives of gas parameters and functions are approximated by the least square method to calculate convective and viscous fluxes in the Navier – Stokes system of equations. The convective fluxes are calculated by the HLLC method. The third-order MUSCL reconstruction scheme is used to achieve high order accuracy. The viscous fluxes are calculated by the second order approximation scheme. The streamlined body surface is represented by a model of an isothermal wall. It implements the conditions for the zero velocity and zero pressure gradient, which is also modeled using the least squares method.

Every moving body is surrounded by its own cloud of points belongs to body’s domain and moving along with it in space. The explicit three-sage Runge–Kutta method is used to solve numerically the system of gas dynamics equations in the main coordinate system and local coordinate systems of each particle.

Two methods for the moving objects modeling with reverse impact on the gas flow have been implemented. The first one uses stationary point clouds with fixed neighbors within the same domain. When regions overlap, some nodes of one domain, for example, the boundary nodes of the particle domain, are excluded from the calculation and filled with the values of gas parameters from the nearest nodes of another domain using the least squares approximation of gradients. The internal nodes of the particle domain are used to reconstruct the gas parameters in the overlapped nodes of the main domain. The second method also uses the exclusion of nodes in overlapping areas, but in this case the nodes of another domain take the place of the excluded neighbors to build a single connected cloud of nodes. At the same time, some of the nodes are moving, and some are stationary. Nodes membership to different domains and their relative speed are taken into account when calculating fluxes.

The results of modeling the motion of a particle in a stationary gas and the flow around a stationary particle by an incoming flow at the same relative velocity show good agreement for both presented methods.

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

The work deals the implementation of a three-dimensional mathematical model of the nonlinear time-varying temperature control and a numerical method of solving it.

Проведено математическое моделирование нестационарных режимов естественной конвекции и поверхностного излучения в замкнутой вращающейся квадратной полости. Рассматриваемая область решения имела две противоположные изотермические стенки, поддерживаемые при постоянных низкой и высокой температурах, остальные стенки являлись адиабатическими. Стенки считались диффузно-серыми. Анализируемая полость вращалась с постоянной угловой скоростью относительно оси, проходящей через центр полости и ориентированной ортогонально области решения. Математическая модель, сформулированная в безразмерных преобразованных переменных «функция тока – завихренность скорости» на основе приближений Буссинеска и диатермичности рабочей среды, была реализована численно методом конечных разностей. Уравнения дисперсии завихренности и энергии решались на основе локально-одномерной схемы А. А. Самарского. Диффузионные слагаемые аппроксимировались центральными разностями, конвективные — с использованием монотонной аппроксимации А. А. Самарского. Разностные уравнения решались методом прогонки. Разностное уравнение Пуассона для функции тока решалось отдельно с применением метода последовательной верхней релаксации. Оптимальное значение параметра релаксации подбиралось на основе вычислительных экспериментов. Анализ радиационного теплообмена проведен с использованием метода сальдо в варианте Поляка. Разработанный вычислительный код был протестирован на множестве сеток, а также верифицирован путем сопоставления полученных результатов при решении модельной задачи с экспериментальными и численными данными других авторов.

Численные исследования нестационарных режимов естественной конвекции и поверхностного теплового излучения в замкнутой вращающейся полости проведены при следующих значениях безразмерных параметров: Ra = 10³–10⁶, Ta = 0–10⁵, Pr = 0.7, ε = 0–0.9. Все распределения были получены для двадцатого полного оборота полости, когда наблюдается установление периодической картины течения и теплопереноса. В результате анализа установлено, что при малой угловой скорости вращения полости возможна интенсификация течения, а дальнейший рост скорости вращения приводит к ослаблению конвективного течения. Радиационное число Нуссельта незначительно изменяется при варьировании числа Тейлора.

Mathematical simulation of unsteady natural convection and thermal surface radiation within a rotating square enclosure was performed. The considered domain of interest had two isothermal opposite walls subjected to constant low and high temperatures, while other walls are adiabatic. The walls were diffuse and gray. The considered cavity rotated with constant angular velocity relative to the axis that was perpendicular to the cavity and crossed the cavity in the center. Mathematical model, formulated in dimensionless transformed variables “stream function – vorticity” using the Boussinesq approximation and diathermic approach for the medium, was performed numerically using the finite difference method. The vorticity dispersion equation and energy equation were solved using locally one-dimensional Samarskii scheme. The diffusive terms were approximated by central differences, while the convective terms were approximated using monotonic Samarskii scheme. The difference equations were solved by the Thomas algorithm. The approximated Poisson equation for the stream function was solved by successive over-relaxation method. Optimal value of the relaxation parameter was found on the basis of computational experiments. Radiative heat transfer was analyzed using the net-radiation method in Poljak approach. The developed computational code was tested using the grid independence analysis and experimental and numerical results for the model problem.

Numerical analysis of unsteady natural convection and thermal surface radiation within the rotating enclosure was performed for the following parameters: Ra = 10³–10⁶, Ta = 0–10⁵, Pr = 0.7, ε = 0–0.9. All distributions were obtained for the twentieth complete revolution when one can find the periodic behavior of flow and heat transfer. As a result we revealed that at low angular velocity the convective flow can intensify but the following growth of angular velocity leads to suppression of the convective flow. The radiative Nusselt number changes weakly with the Taylor number.

В работе развивается иерархический метод математического и компьютерного моделирования интервально-стохастических тепловых процессов в сложных электронных системах различного назначения. Разработанная концепция иерархического структурирования отражает как конструктивную иерархию сложной электронной системы, так и иерархию математических моделей процессов теплообмена. Тепловые процессы, учитывающие разнообразные физические явления в сложных электронных системах, описываются системами стохастических, нестационарных и нелинейных дифференциальных уравнений в частных производных, и в силу этого их компьютерное моделирование наталкивается на значительные вычислительные трудности даже с применением суперкомпьютеров. Иерархический метод позволяет избежать указанных трудностей. Иерархическая структура конструкции электронной системы в общем случае характеризуется пятью уровнями: 1 уровень — активные элементы ЭС (микросхемы, электро-, радиоэлементы); 2 уровень — электронный модуль; 3 уровень — панель, объединяющая множество электронных модулей; 4 уровень — блок панелей; 5 уровень — стойка, установленная в стационарном или подвижном помещении. Иерархия моделей и моделирования стохастических тепловых процессов строится в порядке, обратном иерархической структуре конструкции электронной системы, при этом моделирование интервально-стохастических тепловых процессов осуществляется посредством получения уравнений для статистических мер. Разработанный в статье иерархический метод позволяет учитывать принципиальные особенности тепловых процессов, такие как стохастический характер тепловых, электрических и конструктивных факторов при производстве, сборке и монтаже электронных систем, стохастический разброс условий функционирования и окружающей среды, нелинейные зависимости от температуры факторов теплообмена, нестационарный характер тепловых процессов. Полученные в статье уравнения для статистических мер стохастических тепловых процессов представляют собой систему 14-ти нестационарных нелинейных дифференциальных уравнений первого порядка в обыкновенных производных, решение которых легко реализуется на современных компьютерах существующими численными методами. Рассмотрены результаты применения метода при компьютерном моделировании стохастических тепловых процессов в электронной системе. Иерархический метод применяется на практике при тепловом проектировании реальных электронных систем и создании современных конкурентоспособных устройств.

A hierarchical method of mathematical and computer modeling of interval-stochastic thermal processes in complex electronic systems for various purposes is developed. The developed concept of hierarchical structuring reflects both the constructive hierarchy of a complex electronic system and the hierarchy of mathematical models of heat exchange processes. Thermal processes that take into account various physical phenomena in complex electronic systems are described by systems of stochastic, unsteady, and nonlinear partial differential equations and, therefore, their computer simulation encounters considerable computational difficulties even with the use of supercomputers. The hierarchical method avoids these difficulties. The hierarchical structure of the electronic system design, in general, is characterized by five levels: Level 1 — the active elements of the ES (microcircuits, electro-radio-elements); Level 2 — electronic module; Level 3 — a panel that combines a variety of electronic modules; Level 4 — a block of panels; Level 5 — stand installed in a stationary or mobile room. The hierarchy of models and modeling of stochastic thermal processes is constructed in the reverse order of the hierarchical structure of the electronic system design, while the modeling of interval-stochastic thermal processes is carried out by obtaining equations for statistical measures. The hierarchical method developed in the article allows to take into account the principal features of thermal processes, such as the stochastic nature of thermal, electrical and design factors in the production, assembly and installation of electronic systems, stochastic scatter of operating conditions and the environment, non-linear temperature dependencies of heat exchange factors, unsteady nature of thermal processes. The equations obtained in the article for statistical measures of stochastic thermal processes are a system of 14 non-stationary nonlinear differential equations of the first order in ordinary derivatives, whose solution is easily implemented on modern computers by existing numerical methods. The results of applying the method for computer simulation of stochastic thermal processes in electron systems are considered. The hierarchical method is applied in practice for the thermal design of real electronic systems and the creation of modern competitive devices.

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

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

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

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

При численном решении трехмерных математических моделей ионосферной плазмы используется следующая схема расщепления по физическим процессам: первый шаг расщепления осуществляет продольный перенос, второй шаг расщепления осуществляет поперечный перенос. Исследуемая в работе конечно-разностная схема второго порядка точности приближенно решает уравнения поперечного пере- носа. Эта схема строится с помощью нелинейной процедуры монотонизации Z-схемы, которая является одной из схем второго порядка точности. При этой монотонизации используется нелинейная коррекция по так называемым «косым разностям». «Косые разности» содержат узлы расчетной сетки, относящиеся к разным слоям времени.

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

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

1. Для разрывного начального профиля лучшие свойства показал ограничитель SuperBee.

2. Для непрерывного начального профиля при больших пространственных шагах лучше ограничитель SuperBee, а при малых шагах лучше ограничитель Koren.

3. Для гладкого начального профиля лучшие результаты показал ограничитель Koren.

4. Гладкий ограничитель F показал результаты, аналогичные Koren.

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

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

7. Построенная разностная схема при шагах по времени, не превышающих шаг Куранта, является монотонной и показывает хорошие характеристики точности для решений разных типов. При превышении шага Куранта схема остается устойчивой, но становится непригодной для задач неустойчивости, поскольку условия монотонности перестают в этом случае выполняться.

The purpose of the paper is a research of a 2nd order finite difference scheme based on the Z-scheme. This research is the numerical solution of several two-dimensional differential equations simulated the incompressible medium convection.

One of real tasks for similar equations solution is the numerical simulating of strongly non-stationary midscale processes in the Earth ionosphere. Because convection processes in ionospheric plasma are controlled by magnetic field, the plasma incompressibility condition is supposed across the magnetic field. For the same reason, there can be rather high velocities of heat and mass convection along the magnetic field.

Ionospheric simulation relevant task is the research of plasma instability of various scales which started in polar and equatorial regions first of all. At the same time the mid-scale irregularities having characteristic sizes 1–50 km create conditions for development of the small-scale instabilities. The last lead to the F-spread phenomenon which significantly influences the accuracy of positioning satellite systems work and also other space and ground-based radio-electronic systems.

The difference schemes used for simultaneous simulating of such multi-scale processes must to have high resolution. Besides, these difference schemes must to be high resolution on the one hand and monotonic on the other hand. The fact that instabilities strengthen errors of difference schemes, especially they strengthen errors of dispersion type is the reason of such contradictory requirements. The similar swing of errors usually results to nonphysical results at the numerical solution.

At the numerical solution of three-dimensional mathematical models of ionospheric plasma are used the following scheme of splitting on physical processes: the first step of splitting carries out convection along, the second step of splitting carries out convection across. The 2nd order finite difference scheme investigated in the paper solves approximately convection across equations. This scheme is constructed by a monotonized nonlinear procedure on base of the Z-scheme which is one of 2nd order schemes. At this monotonized procedure a nonlinear correction with so-called “oblique differences” is used. “Oblique differences” contain the grid nodes relating to different layers of time.

The researches were conducted for two cases. In the simulating field components of the convection vector had: 1) the constant sign; 2) the variable sign. Dissipative and dispersive characteristics of the scheme for different types of the limiting functions are in number received.

The results of the numerical experiments allow to draw the following conclusions.

1. For the discontinuous initial profile the best properties were shown by the SuperBee limiter.

2. For the continuous initial profile with the big spatial steps the SuperBee limiter is better, and at the small steps the Koren limiter is better.

3. For the smooth initial profile the best results were shown by the Koren limiter.

4. The smooth F limiter showed the results similar to Koren limiter.

5. Limiters of different type leave dispersive errors, at the same time dependences of dispersive errors on the scheme parameters have big variability and depend on the scheme parameters difficulty.

6. The monotony of the considered differential scheme is in number confirmed in all calculations. The property of variation non-increase for all specified functions limiters is in number confirmed for the onedimensional equation.

7. The constructed differential scheme at the steps on time which are not exceeding the Courant's step is monotonous and shows good exactness characteristics for different types solutions. At excess of the Courant's step the scheme remains steady, but becomes unsuitable for instability problems as monotony conditions not satisfied in this case.