=Paper=
{{Paper
|id=Vol-1482/238
|storemode=property
|title=Матричные методы в задачах механики жидкости
(Matrix methods for fluid mechanics problems)
|pdfUrl=https://ceur-ws.org/Vol-1482/238.pdf
|volume=Vol-1482
}}
==Матричные методы в задачах механики жидкости
(Matrix methods for fluid mechanics problems)==
https://ceur-ws.org/Vol-1482/238.pdf
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
Матричные методы в задачах механики жидкости
М.К. Ермаков
Институт проблем механики им. А.Ю. Ишлинского РАН
Матричный подход используется для моделирования и исследования устойчивости
осесимметричных течений несжимаемой жидкости. Для решения системы нелиней-
ных дискретизованных уравнений Навье-Стокса используется метод Ньютона в мат-
ричном виде. Для исследования устойчивости течения относительно трехмерных
возмущений используется метод обратной итерации и метод Арнольди. Решение сис-
тем линейных уравнений реализуется прямым мультифронтальным методом, реали-
зованным в пакете MUMPS. На основе применения матричного метода для числен-
ного бенчмарка для модели роста кристаллов по Чохральскому удалость получить
решение на сверхподробной сетке с 125M неизвестных при использовании 65 про-
цессоров с памятью 24 ГБ на СК «Ломоносов».
1. Введение
Современную науку и технику невозможно представить без искусственно выращиваемых
кристаллов. Кремний и германий используются в радиоэлектронике, в производстве микро-
схем, создании солнечных батарей. Рубины используются в лазерах и всевозможных датчиках
напряжения, давления, температуры, электромагнитного излучения в различных диапазонах.
Качество выращиваемого из расплава кристалла во многом определяется данным началь-
ным этапом: отсутствием пространственных макронеоднородностей, связанных с неустойчиво-
стью течения в расплаве, и микронеоднородностью распределения примеси, связанной с коле-
баниями фронта кристаллизации [1].
Основой для процессоров является подложка, качество которой обеспечивается устойчи-
вым процессом роста кристалла из расплава. Подавляющее число кристаллов выращивается
методом Чохральского (рис. 1), в котором кристалл вытягивается из расплава, находящимся в
тигле, причем распределение температуры в установке подобрано таким образом, что кристал-
лизация происходит у верхней поверхности расплава. Качество кристалла зависит от кривизны
фронта, однородности распределения легирующей примеси и отсутствия структурных микро-
дефектов. В связи с тем, что легирующая примесь отторгается фронтом, достаточный уровень
конвекции в расплаве необходим для поддержания её равномерной концентрации. Однако,
слишком сильная конвекция приводит к колебаниям и неоднородности.
Рис. 1. Схема гидродинамической модели метода Чохральского роста кристаллов из расплава: 1-
вращающийся диск (кристалл), 2-расплав, 3-тигель.
Целью данной работы является определение границ устойчивости течения расплавов в
представленных конфигурациях роста кристаллов [2-4] на сверхподробных сетках, недоступ-
238
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
ных для расчетов на персональных компьютерах. Моделирование ведется на основе уравнений
Навье-Стокса в приближении Буссинеска. Для решения задачи используются матричные мето-
ды, позволяющие определять основное течение с помощью метода Ньютона и его устойчивость
относительного малых возмущений. Задача об устойчивости течения сводится к решению
обобщенной задачи на собственные значения. Решение возникающих на обоих шагах разреже-
ных линейных систем уравнений осуществляется на основе LU декомпозиции с помощью от-
крытого пакета MUMPS, поддерживающего многопроцессорный режим на основе MPI.
На основе проведенных расчетов для бенчмарка по росту кристаллов по методу Чохраль-
ского для больших чисел Прандтля [3, 4] получены решения на сетках до 5000×5000 точек.
2. Постановка задачи
Гидродинамическая модель роста кристаллов по Чохральскому представляет собой цилин-
дрический объем жидкости (расплава) высоты H и радиуса Rc (рис. 1). Боковая и нижняя сторо-
ны ограничены тиглем, с верхней стороны при r < Rx находится диск (кристалл), вне диска для
Rx < r < Rc поверхность является свободной.
Поведение жидкости определяется уравнениями Навье-Стокса в приближении Буссинеска
для ньютоновской жидкости, поверхностное натяжение которой предполагается линейной
функцией температуры T :
0 (T T0 ) .
Безразмерные уравнения неразрывности, движения и притока тепла в цилиндрических коорди-
натах (r , , z ) с компонентами скорости и давлением P имеют вид:
где
Граничные на твердых поверхностях имеют вид:
Граничные на свободной поверхности имеют вид:
Безразмерные параметры включают два безразмерных геометрических параметра:
, Ar=
а также термокапиллярное число Рейнольдса Re , число Грасгофа Gr ,число Прандтля Pr и чис-
ло Био Bi.
На первом этапе ищутся стационарные осесимметричные решения выписанных уравнений.
На втором этапе найденное решение рассматривается как базовое и для него
рассматривается устойчивость относительно малых возмущений, которые представляются в
виде:
239
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
Относительно амплитудных функций получается обобщенная задача на соб-
ственные значения:
где
Аналогично преобразовываются граничные условия. Описание методики и тесты представлены
в [5].
3. Метод решения и тесты
Для численного решения уравнений используется метод контрольного объема с использо-
ванием разнесенных сеток. При пространственной дискретизации радиальная и осевая компо-
ненты скорости определяются на гранях сеточных объемов, а скалярные величины, давление и
температура, а также азимутальная компонента скорости определяются в центрах этих объе-
мов. Конвективые производные аппроксимируются центральными разностями, обеспечиваю-
щими второй порядок аппроксимации. Введение, при необходимости, существенно неравно-
мерной сетки позволяет воспроизвести тонкую структуру пограничных слоев вблизи точек со-
прикосновения твердых и свободных поверхностей без введения регуляризации.
Для решения стационарных осесимметричных разностных (дискретизованных) уравнений
на первом шаге используется итерационный метод Ньютона в матричной форме для f ( x ) 0 :
f ( x[ k ] ) ,
x[ k 1] x[ k ] k 0,1,...
f ( x[ k ] )
для некоторого начального значения x[0] . Под вектором x подразумевается все множество не-
известных задачи, а под функцией f все множество уравнений и граничных условий. Так как
якобиан f ( x[ k ] ) является матрицей, то для нахождения решения на следующей итерации тре-
буется решать систему линейных алгебраических уравнений.
Использование матричных методов решения уравнений имеет ряд преимуществ по сравне-
нию с традиционными явными и полунеявными методами:
1. возможность построения эффективных полностью неявных схем без расщепления и об-
ременительных ограничений, связанных с шагом интегрирования по времени.
2. отсутствие необходимости разработки метода решения разностных уравнений; разност-
ные уравнения решаются матричными методами.
3. возможность получения решений, неустойчивых для нестационарных уравнений; это не-
обходимо при исследовании линейной устойчивости решений.
240
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
4. Высокая скорость сходимости при использовании метода Ньютона для решения нели-
нейных уравнений; решение уравнений находится за несколько итераций.
Недостатки матричных методов связаны с необходимостью решения систем линейных
уравнений высокого порядка. Матрица якобиана при использовании разностной аппроксима-
ции имеет блочно-диагональный характер и является сильно разреженной. Число элементов
матрицы примерно в 10 раз больше размерности вектора неизвестных. В данной работе для
решения систем линейных уравнений используются прямой метод решения уравнений, осно-
ванный на на реализации мультифронтального метода [6] в пакете MUMPS [7]. Пакет MUMPS
осуществляет LU разложение для разреженных матриц, которое в нашем случае требует при-
мерно в 100 раз больше оперативной памяти, чем исходная матрица. Персонального компьюте-
ра с памятью 12 ГБ достаточно для решения задачи на исходной сетке (nr×nz) до 500×500 то-
чек. Такого разрешения часто вполне достаточно для описания течений вблизи порога их ус-
тойчивости, однако для некоторых случаев этого недостаточно. В частности, для решения бен-
чмарков – специально выделенных тестов для исследования методов и поведения решений.
Решение систем на сетках до 5000×5000 точек в этом случае возможно при использовании кла-
стерной вычислительной системы. При указанном размере сетки количество неизвестных со-
ставит 125M, размер матрицы составит 10 ГБ, а размер LU разложения – 1000 ГБ. Соответст-
венно, задача относится к классу требующей большой оперативной памяти.
Количество даже открытых пакетов для решения задач линейной алгебры весьма велико
[8]. Выбор конкретного пакета следует предпринимать из опубликованного опыта предшест-
венников, что не слишком просто из-за практического отсутствия опыта перекрестного тести-
рования. Вопросы, связанные с выбором программных средств для моделирования задач се-
точными методами на суперкомпьютерах, являются очень важными, и они подробно обсужда-
ются в [9]. Пакет MUMPS является современным структурированным пакетом решения разре-
женных линейных систем, с открытым кодом, поддерживает распределенное хранение матрицы
и LU разложения, основан на шаблонах BLAS и использует MPI [9]. Он также включает стадию
предварительного анализа матрицы для её переупорядочивания с целью уменьшения объема
LU разложения. Помимо встроенных приближенных подходов для переупорядочивания воз-
можно использование специализированных внешних пакетов работы с графами METIS и
SCOTCH. Сочетание данных свойств пакета приводит к увеличению быстродействия до 100 раз
по сравнению с простыми пакетами учебного уровня [10] и уменьшению требований к исполь-
зуемой памяти.
Значительной экономии памяти теоретически можно было бы достичь используя итераци-
онный метод решения градиентного типа для системы линейных уравнений, однако, построе-
ние эффективного предобуславливателя для полной системы уравнений конвекции является
достаточно сложной задачей.
Для решения обобщенной задачи на собственные значения используется метод
обратной итерации и метод Арнольди. Последний используется в режиме shift-and-invert пакета
ARPACK. Для определения устойчивости течения достаточно найти собственное значение с
наибольшей действительной частью. Если действительная части такого собственного значения
положительно, то течение является неустойчивым, в противном случае оно устойчиво. Исполь-
зование стандартных методов решения спектральной задачи затруднено вырожденным харак-
тером матрицы B. Для обоих подходов необходимо снова решение системы линейных уравне-
ний, на этот раз над полем комплексных чисел.
Данный подход был разработан и тщательно тестирован для различных задач вынужден-
ной, термогравитационной и термокапиллярной конвекции [5]. Результаты подтверждены дан-
ными других исследований [11]. Для кластерного варианта кода различие состоит в обращении
к параллельной версии пакета MUMPS (с параллельными пакетами, при необходимости,
ParMETIS и PtSCOTCH для переупорядочивания матриц). Основная программа использует
MPI, синхронно вызывая функции пакета MUMPS и настраивая его на соответствующий режим
работы с помощью входных параметров. Упомянутые пакеты поставляются в исходных текстах
и компиляция библиотек для их использования в среде UNIX не вызывает сложностей. Для
компиляции и выполнения использовались компиляторы фирмы Intel и библиотека MKL.
241
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
4. Результаты
Расчеты проводились на сетках до 5000×5000 точек (125 млн. неизвестных) для нахожде-
ния стационарного решения и до 3000×3000 точек (45 млн. неизвестных) для решения задачи
устойчивости. Для задач использовалось до 65 процессоров с памятью 24 ГБ и до 17 процессо-
ров с памятью 48 ГБ. Результатом решения задачи являются стационарные осесимметричные
течения и нейтральные кривые, определяющие границы устойчивости, критические частоты и
собственные вектора, определяющие наиболее опасные возмущения. (Будут приведены данные
по временным затратам; в данный момент они отсутствуют так расчеты проводились для не-
сколько другой постановки задачи).
В таблице 1 представлено сравнение результатов по бенчмарку для высоких чисел Прандт-
ля [3] для параметров A=0.92, Ar=0.5, Pr=9.2, Re=586, Gr=1.9.105, Bi=0.1. Видна (будет) сходи-
мость результатов при увеличении сеточного разрешения. Результаты, полученные на подроб-
ных сетках могут рассматриваться как практически точные и могут быть использованы для
оценки точности результатов, полученных ранее на менее подробных сетках.
Таблица 1. Сравнение данных по бенчмарку [3]
Настоящая работа [3,4] [3] [4]
Сетка
5000x5000 3000x3000 1000x1000 1000x1000 200x200 120x120
Ψmin -1.5149 -1.5295 -1.510 -1.5175
(r) (0.7959) (0.7969) (0.8100) (0.813)
(z) (0.5468) (0.5557) (0.3263) (0.333)
Ψmax 0.0 0.0 0.0 0.0
Vr min (r=0.5) -41.2904 -52.2848 -36.55 -56.1
(z) (0.9188) (0.9194) (0.9183) (0.9193)
Vr max (r=0.5) 14.6486 14.7508 14.63 14.635
(z) Будет досчитано и пред- (0.08771) (0.08875) (0.08763) (0.08886)
ставлено в окончатель-
Vz min (z=0.5H) -231.33 -232.50 -231.8 -231.5
ном варианте статьи или
(r) ранее по требованию (0.0) (0.0) (0.0) (0.0)
Vz max (z=0.5H) рецензента 15.8016 15.9669 15.76 15.805
(z) (0.9481) (0.9478) (0.9477) (0.9478)
Ekin 81.49 82.66 80.92 81.49
Nu (дно) 1.5860 1.6084 1.575 1.593
Nu (стенка) 1.2559 1.2768 1.246 1.258
Nu (кристалл) -2.8056 -2.8489 -2.784 -2.814
Nu (поверх-
-0.0363 -0.0363 -0.036 -0.0364
ность)
На рис.2 представлены примеры изоповерхностей собственных векторов для бенчмарка [] в
широком диапазоне чисел Прандтля от 0.005 до 30. Отчетливо видно изменение азимутального
волнового числа и формы собственных векторов (возмущений азимутальной компоненты ско-
рости и температуры), определяющей направление распространения возмущений.
242
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
Pr=0.005 Pr=4
Pr=0.05 Pr=10
Pr=0.6 Pr=30
Рис. 2. Примеры собственных функций задачи (температура – слева, азимутальная компонента скорости
- справа) для чисел Прандтля от 0.005 до 30.
5. Выводы
Для осесимметричных конвективных течений роста кристаллов из расплава по методу
Чохральского разработана методика решения стационарных задач и исследования устойчиво-
сти течений по отношению к трехмерным возмущениям для кластерных систем. Для решения
задач используется матричный подход, основанный на применении методов Ньютона, обрат-
ных итераций, Арнольди и линейного решателя (солвера) MUMPS для прямого решения разре-
женных систем.
Использование данного подхода позволило получить численное решение на сетках до
5000×5000 точек (125 млн. неизвестных) на 65 процессорах с памятью 24 ГБ или на 17 процес-
сорах с памятью 48 ГБ. Результаты получены для международного бенчмарка по росту кри-
сталлов из расплава по Чохральскому для высоких чисел Прандтля [3].
Разработанные методы и коды могут быть также использованы для исследования устойчи-
вости течений в задачах внешнего обтекания, крупномасштабной конвекции в мантии и атмо-
сфере Земли и динамики плазмы в осесимметричных конфигурациях.
Работа выполнена с использованием ресурсов суперкомпьютерного комплекса МГУ имени
М.В. Ломоносова [12, 13].
Литература
1. Lappa M. Rotating Thermal Flows in Natural and Industrial Processes. Wiley, Chichester, 2012,
522 p.
2. Wheeler A.A. Four test problems for the numerical simulation of flow in Czochralski crystal
growth // J. Cryst. Growth. 1991. V. 102. № 4. P. 691-695.
243
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
3. Crnogorac N., Wilke H., Cliffe K.A., Gelfgat A.Yu., Kit E. Numerical modelling of instability and
supercritical oscillatory states in a Czochralski model system of oxide melts // Cryst. Res. Tecnol.
2008. V. 43. № 6. P. 606–615.
4. Бессонов О.А., Полежаев В.И. Нестационарные неосесимметричные течения в гидродина-
мической модели метода Чохральского при больших числах Прандтля // Изв. РАН. МЖГ.
201. № 5. С. 16–32.
5. Ermakov M.K., Ermakova M.S. Linear-stability analysis of thermocapillary convection in liquid
bridges with highly deformed free surface // J. Cryst. Growth, 2004. Vol. 266. P.160-166.
6. Amestoy P.R., Duff I.S., L’Excellent J.-Y. Multifrontal parallel distributed symmetric and
unsymmetric solvers. Comput. Methods Appl. Mech. Eng., 184:501–520, 2000.
7. MUMPS. URL: http://mumps-solver.org.
8. Freely Available Software for Linear Algebra.
URL: http://www.netlib.org/utk/people/JackDongarra/la-sw.html.
9. Василевский Ю.В., Коньшин И.Н. Копытов Г.В., Терехов К.М. INMOST – программная
платформа и графическая среда для разработки параллельных численных методов на сетках
общего вида. М.: Изд-во Московского университета, 2013. – 144 с.
10. Saad Y. Iterative methods for sparse linear systems. SIAM. 2003. 528 p.
11. Gelfgat A. Yu. Numerical study of three-dimensional instabilities of Czochralski melt flow driven
by buoyancy convection, thermocapillarity and rotation // Studies on Flow Instabilities in Bulk
Crystal Growth. Transworld Res. Network, 2007. 37/661 (2). P. 1–26.
12. Ермаков М.К. Исследование устойчивости в гидродинамических моделях роста кристаллов
// Cуперкомпьютерные технологии в науке, образовании и промышленности: Альманах / –
М.: Изд-во Московского университета, 2014. С. 140-145.
13. Воеводин Вл.В., Жуматий С.А., Соболев С.И., Антонов А.С., Брызгалов П.А., Никитенко
Д.А., Стефанов К.С., Воеводин Вад.В. Практика суперкомпьютера "Ломоносов" // Откры-
тые системы. – М.: Издательский дом "Открытые системы", N 7, 2012. С. 36-39.
244
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
Matrix methods for fluid mechanics problems
Michael Ermakov
Keywords: supercomputing, matrix methods, Navier-Stokes equations, MUMPS, ARPACK
A matrix approach is used for modeling and stability analysis of axisymmetrical fluid flows
for the Navier-Stokes equations in Boussinesq approximation. Advantages of the matrix
approach consist in an unnecessity of equations splitting and an unnecessity of a partial
method building for equations solution. For the solution of a system of discretized Navier-
Stokes equations the matrix Newton method is employed. At each iterations a linear system
for a vector of the full set of unknowns and for the Jacobian matrix is solved. For a flow
stability analysis against 3-D disturbances, are being considered as a decomposition of normal
modes in azimuthal direction, the inverse iterations and Arnoldi methods (ARPACK) are
used. For a linear-system solution the direct method (MUMPS) is used. Due to application of
matrix approach for numerical benchmark for Czochralski method of crystal growth, it was
possible to get a solution for superfine mesh with 125M unknowns by use of 65 processors
(24 GB) of the “Lomonosov” supercomputer.