Параллельная реализация многомасштабного подхода для расчета течений газа в микроканалах технических систем*

Abstract

Федеральное государственное бюджетное учреждение науки Институт прикладной математики им. М.В. Келдыша Российской академии наук1 Работа посвящена трехмерному моделированию течений газов в микроканалах сложных технических систем. Предложен многомасштабный подход, сочетающий решения уравнений квазигазодинамики (КГД) и молекулярной динамики (МД). Представлена параллельная реализация подхода, основанная на методах расщепления по физическим процессам и разделения областей. Реализация ориентирована на использование вычислительных систем с центральной и гибридной архитектурами. Расчеты показали устойчивость численного алгоритма. С его помощью методами МД были получены основные коэффициентные зависимости для КГД системы, проверен переход от МД к КГД и обратно, произведен расчет трехмерного течения. Полученные результаты подтвердили эффективность разработанного подхода. Ключевые слова: газовая динамика, молекулярная динамика, параллельные алгоритмы и программы, микроканалы, многомасштабные вычисления.

-

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

В рамках данной работы предлагается следующий комбинированный подход. Он сочетает расчеты течения газовой среды в конкретных микроканалах и расчеты взаимодействия струй с подложкой, а также перемещения наночастиц в потоке газа и взаимодействие их с подложкой при попадании к месту напыления. При этом моделирование проводится на двух масштабных уровнях – основном, имеющем характерные размеры от нескольких десятков микрон до десятков миллиметров, и дополнительном, имеющем характерные размеры порядка микрона и менее. На основном уровне (в макромасштабах) рассчитываются течения газа и наночастиц во всех частях установки, а также учитываются внешние факторы, связанные с управлением процессом напыления. На дополнительном уровне (в микромасштабах) рассчитываются взаимодействия: 1) молекул газа между собой (формирующие уравнение состояния смеси и реализующие перемешивание компонент), 2) молекул газа и атомов твердых поверхностей (описывающие явления в погранслоях), 3) молекул газа и атомов наночастиц (описывающие образование погранслоев на поверхности наночастиц), 4) молекул газа и атомов подложки (учитывающие процессы теплообмена с подложкой), 5) атомов наночастиц с атомами подложки (процессы отталкивания, прилипания, внедрения наночастиц в подложку).

Указанные физические процессы являются сложными и многостадийными. Поэтому на данном этапе исследований используется упрощенное описание проблемы. Оно предполагает расщепление по физическим процессам. Фактически основной цикл моделирования представляет собой чередование во времени расчетов на макро- и микроуровне. В представляемом варианте исследования рассматривается расчет течения азот-водородной смеси после выхода из сопла в полуоткрытый микроканал и далее в свободное пространство (см. Рис. 1). Рис. 1. Геометрия модельной задачи, для удобства представленная в виде двумерного среза. Из баллона через микросопло подается в микроканал газовая смесь, использующаяся для управления траекториями наночастиц, которые доставляются в нужные места подложки Некоторая часть исследования была проделана в предыдущих работах [ 2-6 ]. В частности, в [ 2-5 ] развивались методы расчета термодинамического равновесия в системах металл-металл, газ-газ и газ-металл. В работе [ 6 ] было предложено и частично апробировано объединение микро- и макромоделей газовой среды в единый объект исследования.

Основной задачей в данном исследовании была разработка полного трехмерного многомасштабного численного метода и параллельного алгоритма для его реализации пригодного для расчета течения газовой смеси в микроканале и вблизи подложки с учетом граничных эффектов, разрешаемых методами молекулярной динамики [ 7 ]. Описание деталей модели и алгоритма дается в последующих пунктах. 2. Теоретическая часть 2.1 Макроскопическая модель

Макроскопическая модель течения во всех частях установки основывается на квазигазодинамических (КГД) уравнениях [ 8 ]. Выбор КГД уравнений связан с двумя факторами. Вопервых, КГД уравнения хорошо зарекомендовали себя при расчетах сверхзвуковых течений, разреженных вязких и теплопроводных газов. Во-вторых, они хорошо подходят к расчетам течений в микросистемах, так как естественным параметром обезразмеривания в них является средняя длина свободного пробега. В-третьих, КГД система уравнений применима в широком диапазоне чисел Кнудсена, что дает возможность моделировать сложную систему с микроканалами разных диаметров, используя одну и ту же математическую модель.

В случае смеси газов система КГД уравнений записывается для каждого газа отдельно и имеет одинаковый вид [ 8-9 ]. Запишем эти уравнения в случае бинарной смеси в инвариантном относительно системы координат виде вместе с уравнениями связей и состояния: l  div Wl   0, t

 1  Wl   lul  lwl , wl  ul ,ul   pl  , l  a,b;

 l   t lul,k  div Wluk   Sluk  ,

 2  Wluk   lulul,k  ek  pl  3 l divul   l ul,k  ,ek ul   l wl,kul  l wlul,k ,  Sluk   ll'l ul*,k  ul,k , l  a,b, l '  b, a, k  x, y, z;  t El  div WlE  SlE ,

 2  WlE  lul  lwl  Hl  lTl   3 l divul ul   l ul,k  ,ek ul   l wl ,ul ul , kx,y,z SlE  ll'l  El*  El , l  a,b, l '  b,a, (1) (2) (3) El  12 l ul 2  l l , Hl  El  pl , pl  ZlllTl ,  l  сV ,lTl , l  a,b. (4) l Здесь предполагается, что смесь состоит из газов а и b с числовыми плотностями (концентрациями) na и nb и массовыми плотностями  a  mana и  b  mbnb , где ma и mb – массы молекул газов а и b. Каждый газ характеризуется своей температурой Tl и макроскопической скоростью ul . Другие параметры компонент смеси: pl – парциальные давления, El – плотности полной энергии,  l – внутренние энергии, Hl – полные энтальпии, l  l Tl  и l  l Tl  – коэффициенты вязкости и теплопроводности, Zl  Zl Tl  – коэффициенты сжимаемости, сV ,l  сV ,l Tl  – удельные теплоемкости при постоянном объеме, l  kB / ml – газовые постоянные ( kB – постоянная Больцмана). Векторы Wl  , Wluk  , WlE с точностью до знака совпадают с потоками плотности, соответствующих компонент плотности импульса и плотности энергии, l wl – векторы КГД поправки к потоку плотности, пропорциональные максвелловскому времени релаксации  для смеси газов, ek – единичные орты.

Обменные члены Sluk  и SlE учитывают перераспределения импульса и энергии между компонентами смеси. Они содержат компоненты скоростей u*l и энергии El* , вычисляемые на молекулярном уровне (см. п. 2.2), а также  ll ' – частоты взаимных столкновений между молекулами сорта l c молекулами сорта l', которые вычисляются через  l – частоты взаимных столкновений между молекулами одного сорта:

1   a 1  Gab b M a   b 1  Gba a Mb  ,

 a Mb   b M a  Gab  1  a Mb 2  2 21  M a  1 .

 b M a    Mb  

1  dab   ab  a   da  2 ma  mb nb ,  ba  ab nna  b  dab  2mb na b  db  2 ma  mb na , 2ma nb dab  da 2 db ,  l  pll l , l  a,b. где dl – диаметры молекул сорта l, l – безразмерные параметры, характеризующие процесс столкновения молекул сорта l [ 10 ].

Коэффициенты вязкости l , теплопроводности l и сжимаемости Zl , а также удельные теплоемкости сV ,l и средние длины свободного пробега l для компонент смеси определяются либо из базы данных молекулярных расчетов, либо с помощью прямых МД вычислений (см. п. 2.2). Локальные скорости звука для компонент смеси рассчитываются по формулам:  al   l ZllTl ,  l  1 l , l  a,b,

cv,l где  l   l Tl  – показатели адиабат.

Параметры газовой смеси определяются следующим образом:   a  b , u  aua  bub , p  pa  pb , E  Ea  Eb , n  na  nb ,

 T  naZaTa  nbZbTb , m  mana  mbnb , Z  Za  Zb ,   nZ n p  ZT

k   B ,   .

m p Для определения вязкости смеси используем формулу Будденберга-Уилки [ 11 ]: (5) (6) (7) (8) (9) (10) Здесь Ma и Mb – молярные массы газов а и b. Величина Gba вычисляется аналогично путем циклической замены индексов. Средние локальная скорость звука и длина свободного пробега, числа Маха и Рейнольдса для смеси равны: a  maaa  mbab  ma  mb ,   maa  mbb  ma  mb , Ma  u a , Re   u   . Система уравнений (1)-(3) замыкается начальными и граничными условиями. Начальные условия соответствуют равновесному состоянию газовой среды в отсутствие взаимодействия с внешними факторами. Рассмотрен случай покоящейся среды во всех точках расчетной области: l  l,0, ul  0, pl  pl,0, Tl  T0, l  a,b. Здесь l,0 , pl,0 – начальные плотности и давления компонент газа, T0 – начальная температура.

Область сопла и баллон со смесью газов не входили в расчет. Полагалось, что в начальный момент времени открывается перегородка между соплом и областью микроканала, и газ из области высокого давления начинает поступать в область низкого давления. Сопло считалось идеальным адиабатическим соплом Лаваля. Параметры газовой среды на выходе из сопла (то есть на входе в среду) рассчитывались в соответствии с известными формулами:

Здесь параметры с индексом bal соответствуют значениям параметров в баллоне, а параметры с индексом in – соответствующим параметрам на входе в среду.

На стенках микроканала и на подложке ставятся следующие граничные условия:  Wl( ) ,n   l l  l(w) , Wl(uk ) ,n  l,k lul,k  l,wul(,wk) , k  x, y, z,  Wl(E) ,n  l  El  El(w) , l  a,b.

В общем случае молекулы газа могут проникать в металлическую стенку, либо отражаться от нее под некоторым углом. Взаимодействие со стенкой может быть упругим или неупругим, термическим или изотермическим. Удобно вблизи стенки ввести специальную микросистему, состоящую из молекул газа и атомов металла. В этой микросистеме газовые компоненты имеют плотности массы l(w) , импульсов l(w)ul(w) и полной энергии El(w) . Граничные условия на стенке можно записать в виде условий третьего рода, описывающих обмен массой, компонентами импульса и энергией между газовой смесью в потоке и вблизи стенок. Коэффициенты  l , l,k и  l , связанные с коэффициентами аккомодации [ 12 ], априори неизвестны, но могут быть определены либо из таблиц физических величин (что возможно лишь для ограниченного диапазона температур и давлений), либо рассчитаны методами молекулярной динамики. Второй подход представляется наиболее оправданным, поскольку он является универсальным и поможет согласовывать процессы взаимодействия на микро- и макроуровнях.

На свободных поверхностях задаются так называемые "мягкие" граничные условия [ 8 ]: (11) (12) (13) (14) l  0 , n  lul   0 , pl  0, l  a,b.

n n 2.2 Микроскопическая модель

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

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

В данной работе остановимся на случае чистой бинарной смеси и одного сорта металлических частиц. Тогда на микроуровне система газ-металл представляется совокупностью частиц, поведение которых описывается уравнениями Ньютона [ 7 ]:

dvl,i  Fl,i , vl,i  drl,i , i  1,..., Nl , l  a,b,c, ml,i dt dt где i – номер частицы, l  a,b,c – тип частиц ( a – молекулы первого газа, b – молекулы второго газа, c – атомы металла), Nl – полное число частиц типа l , ml,i – масса частицы типа l с номером i , rl,i  rx,l,i , ry,l,i , rz,l,i  и vl,i  vx,l,i ,vy,l,i ,vz,l,i  – радиус-вектор и вектор скорости i -ой частицы типа l , Fl,i   Fx,l,i , Fy,l,i , Fz,l,i  – суммарная сила, действующая на данную частицу.

Силы представляют собой сумму компоненты взаимодействия i -ой частицы с окружающими частицами и компоненты, отвечающей за внешнее воздействие:

Fl,i   U rl,1,...,rl,Nl   Fle,xit , i  1,..., Nl , l  a,b,c,

rl,i где U – суммарная потенциальная энергия, Fle,xit – сила взаимодействия с внешней средой.

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

U  Uaa  Ubb  Ucc  Uab  Uac  Ubc ,

Na Nb Uaa   aa  ra,i  ra, j , Ubb   bb  rb,i  rb, j , i1, ji i1, ji

Nc   Ucc   1,cc rc,i    2,cc  rc,i  rc, j  , Uab  i1  ji  1 Na Nb

  ab  ra,i  rb, j , 2 i1 j1 Uac  1 Na Nc

  ac  ra,i  rc, j , Ubc  2 i1 j1 1 Nb Nc

 bc  rb,i  rc, j .

2 i1 j1 Каждый вид взаимодействия описывается с помощью соответствующего потенциала ll' . Для смеси водород-азот использовался потенциал Ми в форме "n-6" [ 13 ], адаптированный в работе [ 14 ] к расчетам смеси водорода и азота в зависимости от их объемного соотношения в выделенном микрообъеме. Для взаимодействия атомов никеля между собой использовалась форма потенциала EAM [ 15 ], который учитывает не только парные взаимодействия 2,cc , но и влияние окружения на определенную частицу 1,cc . Для учета взаимодействий газ-металл использовались стандартный потенциал Леннард-Джонса [ 16 ] и потенциал Морзе [ 17 ].

Вычисление коэффициентов сжимаемости и теплоемкости подробно описано в работе [ 4 ]. Коэффициенты вязкости и теплопроводности вычислялись по формулам [ 18 ]: l  1 

  p,l t0  t   p,l t0  dt ,   xy, yz, zx, l  a,b, kBTlV 0 p,l   ml,i  v,l,i  v,l,i   r,l,i  r,l, j    F,l,i  F,l, j , i, j  1...Nl , i i ji 1  l  kBTl2V 0 q,l t0  t   q,l t0  dt ,   x, y, z, q,l  ddt,l , ,l  i r,l,i l,i  l,i , l,i  ml 2vl,i2  1 N

ll  rl,i  rl, j , 2 ji где V – объем исследуемой системы, l,i – мгновенная полная энергия молекулы i типа l , t0 – начальный момент времени при расчетах коэффициентов.

Средние длины свободного пробега вычисляются по формулам [ 19 ]: (15) (16) (17)  m   m  a  1  nada2 2  nbda2b 1  a , b  1  nbdb2 2  nada2b 1  b ,  mb   ma  (18) где l – эффективный радиус взаимодействия молекул сорта l [ 14 ].

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

Граничные условия на молекулярном уровне выбираются в зависимости от моделируемой ситуации. В частности, когда расчеты направлены на определение свойств среды, рассматривается некоторый выделенный трехмерный ее объем, вне которого предполагается периодическое продолжение среды на неограниченные расстояния по всем трем пространственным координатам. В этом случае используются периодические граничные условия [ 7 ].

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

В данной работе использовались и периодические, и зеркальные граничные условия, а также вход потока частиц в среду и свободный выход частиц из нее. Дополнительно отдельные части микросистем термостатировались [ 20, 21 ]. 3. Реализационная часть

Реализация разработанного многомасштабного подхода опиралась на численный алгоритм, основанный на расщеплении по физическим процессам и использующий сеточные аппроксимации КГД уравнений (1)-(3) и подсеточные МД вычисления по уравнениям (14)-(16).

Алгоритм состоит в следующем. В начальный момент времени на выбранной сетке задается равновесное состояние макросистемы (10). При необходимости вблизи границ расчетной области задается также равновесное состояние газовой и металлической микросистем.

Далее на каждом шаге по времени производится сначала предикторный расчет макропараметров компонент газа по сеточным аналогам КГД уравнений (1)-(3) без учета обменных членов (КГД блок программы). В результате расчета определяются новые макропараметры компонент газа в каждом контрольном объеме пространственной сетки.

Затем производится подсеточный МД расчет с целью вычисления обменных членов и кинетических коэффициентов (МД блок программы). Он осуществляется независимо в каждом контрольном объеме (ячейке) сетки и производится с существенно более мелким шагом по времени, связанным с эволюцией молекулярной подсистемы. Критерием останова МД расчета является либо достижение характерного времени эволюции молекулярной системы, либо изменение (на 1-2%) одного или нескольких макропараметров молекулярной системы (средние импульс, кинетическая или потенциальная энергии). Если сильных изменений макропараметров не происходит, то расчет производится до достижения некоторого заданного момента времени пропорционального времени максвеллизации МД системы. Там же в МД блоке осуществляется возврат к макроскопическому уровню и проводится коррекция плотностей импульсов и плотностей полной энергии с помощью вычисленных обменных членов. Коррекция проводится в каждой точке сетки. Завершается процесс МД вычислений расчетом кинетических коэффициентов и уравнений состояния. Таким образом, связь макро- и микроуровней в сеточной пространственной и временной областях осуществляется посредством поочередного пересчета макропараметров по КГД и МД уравнениям. В случае присутствия в потоке наночастиц используется третий блок вычислений (НЧ блок программы), в котором решаются либо уравнения конвекции-диффузии для концентраций наночастиц каждого сорта (эта ситуация рассматривается для случая мелких наночастиц), либо уравнения Ньютона для баллистического переноса наночастиц (в случае крупных наночастиц). В данной работе эта часть вычислений не рассматривается.

Завершается описанная часть алгоритма проверкой критерия окончания расчетов и переходом к следующему шагу по времени в случае невыполнения критерия. 3.1 Конечно-объемная схема решения КГД уравнений

Для расчета макропараметров по КГД уравнениям разработан сеточный численный метод, имеющий в своей основе метод контрольных объемов [ 22-24 ]. Для этого в расчетной области D была введена пространственная сетка D с ячейками Cm ( m 1,..., M ) и сетка по времени  t с переменным шагом t . Сетка D в общем случае является гибридной, то есть включает в себя несколько типов элементов: тетраэдры, пентаэдры, гексаэдры, гептаэдры и октаэдры.

Все параметры компонент газа (плотности, давления, температуры, компоненты вектора скорости и т.д.) были отнесены к центрам масс элементов сетки, то есть к центрам ячеек. Потоковые переменные были заданы в центрах граней ячеек. Пространственные аппроксимации основных членов были выполнены по методике, представленной в работе [ 25 ]. Вычислительная схема по времени была явной. Итоговые сеточные уравнения на этапе предиктора имеют вид: ˆl,m  l,m  t Mm

Wl(,n,)m, p Sm, p , Wl(,n,)m, p  Wl(,m),p ,nm,p , Vm p1 t Mm

Wl(,n,um), p,k Sm, p , Wl,(n,um), p,k   Wl(,mu,kp) ,nm,p , k  x, y, z, Vm p1 (19) ˆl,muˆl,k,m  l,mul,k,m  ˆ El,m  El,m  t Mm

Wl(,nE,m), p Sm, p , Wl(,nE,m), p  Wl(,Em), p ,nm, p , m  1,..., M , l  a,b.

Vm p1 Здесь fm , fˆm – значения соответствующих функций на нижнем и верхнем слоях по времени в центре контрольного объема Vm (совпадающего в силу выбранной аппроксимации с ячейкой сетки Cm ), Vm – величины контрольных объемов Vm , Sm, p – площади граней Sm, p контрольного объема Vm , M m – количество граней контрольного объема Vm , Wn(,m), p – скалярные произведения потока

Wm(, p) через грань Sm, p на внешнюю нормаль nm, p к этой грани (   ,ux ,uy ,uz , E ). Аппроксимации потоков Wm(, p) аналогичны предложенным в [ 25 ], но не используют искусственную вязкость.

Полученные предикторные переменные ˆl , ˆluˆ l , Eˆl передаются в блок молекулярной динамики. В результате его работы вычисляются модифицированные параметры ˆluˆ *l и Eˆl* . По ним находятся окончательные значения плотностей импульсов ˆluˆ *l* и энергии Eˆl** : ˆl,muˆl*,*k,m  ˆluˆl,k,m  t ll'ˆl,m uˆl,k,m  uˆl*,k,m , k  x, y, z,

* Eˆl*,*m  Eˆl*,m  t ll'ˆl,m  Eˆl,m  Eˆl*,m , m  1,..., M , l  a,b, l '  b,a. (20) По полученным макропараметрам в МД блоке вычисляются все кинетические коэффициенты ( l , l ), длины свободного пробега l и параметры в уравнениях состояния Zl , сV ,l . 3.2 Алгоритм решения уравнений Ньютона Система уравнений (15)-(17) решается с помощью схемы Верле [ 26 ] в скоростной форме: 0) Fl0,i  F ra0,1,...,ra0,Na ;rb0,1,...,rb0,Nb ;rc0,1,...,rc0,Nc , i  1,..., Nl , l  a,b, c; n  0,1, 2,...,:

n 1) rln,i1  rln,i  vln,it  Fl,i (t)2 , i  1,..., Nl , l  a,b,c;

ml 2 3) vln,i1  vln,i  Fln,i1  Fl,i t, i  1,..., Nl , l  a,b,c.

n 2 ml (21) Здесь t – шаг интегрирования, n – номер шага, Fln,i – значение силы на шаге n , F – процедура вычисления сил на основе формул (15), (16).

Для достижения желаемых температур газа и металла используется термостат Берендсена [ 20 ], для достижения желаемых температуры и импульса газа – термостат Ланжевена [ 21 ]. 3.3 Параллельная реализация

Параллельная реализация алгоритма предполагает использование кластера (или суперкомпьютера) с центральной или гибридной архитектурой, имеющего на каждом узле несколько многоядерных центральных процессоров (ЦПУ), а также несколько векторных или графических процессоров (ВПУ или ГПУ). Распараллеливание алгоритма производится на принципах геометрического параллелизма и разделения областей. Основной газодинамический расчет производится по дискретным КГД уравнениям на сетке, распределенной между узлами кластера с помощью техники "domain decomposition" [ 27 ]. Внутри узла КГД вычисления распределяются между потоками ЦПУ. Подсеточные МД вычисления возлагаются на ВПУ или ГПУ при их наличии. Распараллеливание КГД вычислений между потоками ЦПУ также производится геометрическим способом. Распараллеливание МД вычислений производится путем разбиения всего множества частиц, относящихся к одной ячейке сетки, на группы одинаковой мощности (подробно эти расчеты представлялись в [ 5 ]). Каждый блок трэдов ВПУ или ГПУ обрабатывает одну или несколько молекулярных групп, относящихся к одному или нескольким узлам сетки. Программная реализация алгоритма выполнена по гибридной технологии [ 28 ], использующей MPI [ 29 ], OpenMP [ 30 ] и CUDA [ 31 ]. 4. Экспериментальная часть

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

В указанной расчетной области была построена цилиндрически симметричная призматическая сетка, полученная прямым произведением двумерной треугольной сетки в поперечном сечении канала и равномерной сетки по продольной координате. Прямые треугольные призмы имели в основании правильный треугольник со стороной a(k)  a(0) / 2k , и высоту также равную a(k) ( a(0) =1.5 мкм, k  0,1,2,3,... – параметр измельчения сетки). Стартовая сетка (0) имела M (0)  6 25 60  9000 ячеек. Последующие сетки получались измельчением ячеек предыду

V щей сетки на 8 частей. В итоге сетка (k) имела MV(k)  MV(0)8k ячеек. Объем каждой ячейки в такой сетке равен V (k)  0.5 3 a(k) 3 мкм3. Количество ячеек, примыкающих к боковой поm каждой ячейке Cm(0) стартовой сетки (0) (имеющей объем V (0)  2.3864851018 м3) при н.у. m ( p 101325 Па , T  273.15К ) содержится N p(0a)rt  NA Vm(0) V1  6.412107 молекул газа ( N A – число Авогадро, V – объем моля газа при н.у.). Этого количества достаточно, чтобы получить представительную статистику при МД вычислениях, однако, слишком много для реальных вычислений. Поэтому в длительных расчетах лучше использовать более подробные сетки, содержащие большее число ячеек, но имеющие меньший объем одной ячейки и содержащие меньшее число частиц в ней. Как показал анализ, оптимальное соотношение количества ячеек сетки и частиц в одной ячейке сетки достигается для сетки (3) , содержащей 4608000 ячеек, в каждой из которых при н.у. содержится 125234 молекул газа. Если начальное давление в канале составляет 0.001 атм, то стартовое количество молекул газа в ячейке падает до 125.

Если металлическая поверхность рассматривается на атомном уровне, то в вычислениях МД блока появляется необходимость учета структуры кристаллической решетки никеля. Как было показано в работе [ 3 ], при н.у. размер ребра гранецентрированной решетки никеля равен aNi  0.35314 нм . При толщине учитываемого металлического слоя в 4 кристалла, в ячейках, примыкающих к боковой поверхности канала, необходимо учесть взаимодействие молекул газа 2 с 4a(k) / aNi  кристаллами; в ячейках, примыкающих к левому торцу канала, – взаимодействие с примерно вдвое меньшим числом кристаллов. Для сетки (0) в этом случае получаем соответственно 72168667 и 36084333 кристаллов, что потребует неоправданно больших вычислительных затрат. Для оптимальной сетки (3) эти величины составляют 1127635 и 563817, что более реалистично. Однако еще лучше в этом случае использовать сетку (4) .

Тестирование разработанного подхода проводилось на суперкомпьютерах К-100 (ИПМ им. М.В. Келдыша РАН), МВС-10П (МСЦ РАН), кластере с сетью Ангара-К1 (ОАО "НИЦЭВТ") и кластере IBM Power8 (ЛИТ ОИЯИ). Расчеты проводились на центральных (ЦПУ – Intel Xeon X5670, Intel Xeon E5-2690, Intel Xeon E5-2630 и IBM Power E870), векторных (ВПУ – Intel Xeon Phi 5110p) и графических (ГПУ – NVidia K40m) процессорах. Количество ядер ЦПУ на разных системах варьировалось от 16 до 2048. Количество ядер ВПУ варьировалось от 60 до 7680. Количество видеоядер ГПУ варьировалось от 2880 до 11520. Расчеты проводились как без учета атомного строения металлической поверхности, так и с учетом этого фактора. В первом случае использовалась сетка (3) , во втором – (4) .

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

Во-вторых, общее время расчета трехмерной задачи существенным образом зависит от выбранной сетки и от соотношения времен расчетов КГД и МД блоков. Если сетка крупная, то основная вычислительная емкость сосредоточена в МД блоке, и расчеты в этом случае существенно замедляются. Если сетка достаточно подробная, то емкости КГД и МД блоков могут быть сравнимы. Однако при этом следует обращать внимание на количество молекул газа в ячейках, чтобы иметь представительную статистику для получения макропараметров газа.

В-третьих, эффективность распараллеливания в расчетах на ЦПУ оказалось в 1.5 раза выше, чем в расчетах на ВПУ. Максимальная эффективность распараллеливания на 2048 ядрах ЦПУ достигла 83.2% при решении тестовой задачи на сетке (4) . Для расчетов на сетке (3) она была порядка 64.8%, однако здесь следует учесть, что вычисления во втором случае не учитывают строение металлической поверхности канала и выполняются на порядок быстрее.

В-четвертых, при расчетах на ЦПУ и ГПУ скорость расчетов МД блока повысилась почти в 10 раз. Общая эффективность распараллеливания при этом снизилась на 5.8 % ввиду необходимости перекачки данных из памяти ЦПУ в ГПУ и обратно. Следует отметить, что использовалась сетка (3) и в расчетах участвовали 4 ЦПУ и 4 ГПУ. В дальнейшем планируются расчеты на большем числе узлов, содержащих ГПУ. При использовании ГПУ можно проводить расчеты на более крупных сетках, например, на (2) , так как ГПУ производит вычисления почти на порядок быстрее. Единственным ограничением является размер оперативной памяти ГПУ. 5. Заключение

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

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

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

Тестирование разработанного подхода проводилось на суперкомпьютерах К-100 (ИПМ им. М.В. Келдыша РАН), МВС-10П (МСЦ РАН), на кластере Ангара-К1 (ОАО "НИЦЭВТ") и кластере IBM Power8 (ЛИТ ОИЯИ). Расчеты показали, что общий численный алгоритм устойчив к использованию корректирующих течение данных, полученных в результате МД вычислений. С его помощью методами МД были получены основные коэффициентные зависимости для КГД системы, проверен переход от МД к КГД и обратно, произведен расчет трехмерного течения. Тесты на эффективность распараллеливания показали хорошие результаты. Литература Parallel realization of multiscale approach

for calculating the gas flows in microchannels of technical systems*

V.O. Podryga1, S.V. Polyakov1 Keldysh Institute of Applied Mathematics (Russian Academy of Sciences)1 The work is devoted to three-dimensional modeling of gas flows in microchannels of complex technical systems. Multiscale approach combining solutions of quasigasdynamic equations (QGD) and molecular dynamics (MD) is proposed. Parallel implementation of an approach based on the method of splitting into physical processes and the separation of areas is presented. Implementation is oriented for using computer systems with central and hybrid architectures. Calculations have shown the stability of the numerical algorithm. Basic coefficient dependences for the QGD system were obtained by MD method, the transition from the MD to the QGD and back has been tested, three-dimensional flow calculation was made. The results confirmed the efficiency of the developed approach.

References

1. Alkhimov A.P. , Klinkov S.V. , Kosarev V.F. , Fomin V.M. Kholodnoe gazodinamicheskoe napylenie. Teoriya i praktika [Cold gasdynamic spray . Theory and practice] . Moscow: Fizmatlit, 2010 . 536 p.

2. Podryga V.O. Molecular Dynamics Method for Simulation of Thermodynamic Equilibrium // Mathematical Models and Computer Simulations . 2011 . Vol. 3 , No. 3. P. 381 - 388 .

3. Podryga V.O. , Polyakov S.V. Molecular Dynamics Simulation of Thermodynamic Equilibrium Establishment in Nickel // Mathematical Models and Computer Simulations . 2015 . Vol. 7 , No. 5. P. 456 - 466 .

4. Podryga V.O. Opredelenie makroparametrov real'nogo gaza metodami molekulyarnoy dinamiki [Determination of real gas macroparameters by molecular dynamics] // Matematicheskoe modelirovanie [Mathematical Simulation]. 2015 . Vol. 27 , No. 7. P. 80 - 90 .

5. Podryga V.O. , Polyakov S.V. , Puzyrkov D.V. Superkomp 'yuternoe molekulyarnoe modelirovanie termodinamicheskogo ravnovesiya v mikrosistemakh gaz-metall [Supercomputer molecular modeling of thermodynamic equilibrium in gas-metal microsystems] // Vychislitel'nye metody i programmirovanie [Numerical methods and programming]. 2015 . Vol. 16 , No. 1. P. 123 - 138 .

6. Kudryashova T.A. , Podryga V.O. , Polyakov S .V. Modelirovanie techeniy gazovykh smesey v mikrokanalakh [Simulation of gas mixture flows in microchannels] // Vestnik RUDN . Seriya: Matematika. Informatika. Fizika [Bulletin of People's Friendship University of Russia. Series: Mathematics. Computer science. Physics]. 2014 . No. 3. P. 154 - 163 .

7. Frenkel D. , Smit B. Understanding Molecular Simulation . From Algorithm to Applications. NewYork: Academic Press, 2002 . 638 p.

8. Elizarova T.G. Quasi-gas dynamic equations . Springer, 2009 . 286 p.

9. Elizarova T.G. , Zlotnik A.A. , Chetverushkin B.N. On quasi-gasdynamic and quasi-hydrodynamic equations for binary mixtures of gases // Doklady Mathematics. 2014 . Vol. 90 , No. 3. P. 1- 5 .

10. Bird G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows . Oxford, Clarendon Press, 1994 . 478 p.

11. Buddenberg J.W. , and Wilke C.R. Calculation of Gas Mixture Viscosities // Ind. Eng. Chem . 1949 . Vol. 41 , No. 7. P. 1345 - 1347 .

* This work was supported by the Russian Foundation for Basic Researches (projects № 15 - 07 -06082-а, № 15 - 29 -07090-ofi_m)

12. Abramovich G.N. Prikladnaya gazovaya dinamika . Ch. 1 [Applied gas dynamics. Part 1 ]. Moscow: Publishing of Science, 1991 . 600 p.

13. Von Mie G. Zur kinetischen theorie der einatomigen korper // Ann. Phys., Leipzig. 1903 . Vol. 11 , No. 8. P. 657 - 697 .

14. Fokin L.R. , Kalashnikov A.N. The transport properties of an N2-H2 mixture of rarefied gases in the EPIDIF database // High Temperature . 2009 . Vol. 47 , No. 5. P. 643 - 655 .

15. Zhou X.W. , Johnson R.A. , Wadley H.N.G. Misfit-energy-increasing dislocations in vapordeposited CoFe /NiFe multilayers // Phys. Rev. B. 69 . 2004 . P. 144113 .

16. Lennard-Jones J .E. Cohesion // Proceedings of the Physical Society . 1931 . Vol. 43 , No. 5. P. 461 - 482 .

17. Morse P.M. Diatomic molecules according to the wave mechanics . II. Vibrational levels // Phys. Rev . 1929 . Vol. 34 . P. 57 - 64 .

18. Maruyama S. Molecular dynamics method for microscale heat transfer / / In: Advances in Numerical Heat Transfer (eds. W. J. Minkowycz , E. M. Sparrow ) . Taylor & Francis , New York, 2000 . Vol. 2 , Chap . 6. P. 189 - 226 .

19. Sivukhin D . V. Obshchiy kurs fiziki . T. II. Termodinamika i molekulyarnaya fizika [The general course of physics. Vol. II. Thermodynamics and molecular physics] . Moscow: Publishing of Science, 1975 . 519 p.

20. Berendsen H.J.C. , Postma J.P.M. , W.F. van Gunsteren et al. Molecular dynamics with coupling to an external bath // J. Chem. Phys. 1984 . Vol. 81 . P. 3684 - 3690 .

21. Heerman D .W. Computer simulations methods in theoretical physics . Springer, 1986 . 148 p.

22. Fryazinov I.V. Metod balansa i variatsionno-raznostnye skhemy [Balance method and variationaldifference schemes] // Differentsial'nye uravneniya [Differential Equations]. 1980 . Vol. 16 , No. 7. P. 1332 - 1343 .

23. Samarskiy A.A. , Koldoba A.V. , Poveshchenko Yu .A., Tishkin V.F. , Favorskiy A.P. Raznostnye skhemy na neregulyarnykh setkakh [Difference schemes on irregular grids] . Minsk: Publishing of Closed Joint Stock Company “Criterion”, 1996 . 196 p.

24. Eymard R. , Gallouet T.R. , Herbin R . The finite volume method // In: Handbook of Numerical Analysis (eds. P.G. Ciarlet, J.L . Lions). Amsterdam: North Holland, 2000 . Vol. VII. P. 713 - 1020 .

25. Popov I.V. , Fryazinov I.V. Method of adaptive artificial viscosity for the equations of gas dynamics on triangular and tetrahedral grids // Mathematical Models and Computer Simulations . 2013 . Vol. 5 , No. 1. P. 50 - 62 .

26. Verlet L. Computer «experiments» on classical fluids. I. Thermodynamical properties of LennardJones molecules // Phys. Rev . 1967 . Vol. 159 . P. 98 - 103 .

27. Toselli A. , Widlund O. Domain Decomposition Methods - Algorithms and Theory. Springer Series in Computational Mathematics. Springer, Berlin, 2005 . Vol. 34 . 450 p.

28. Polyakov S.V. , Karamzin Yu .N., Kosolapov O.A. , Kudryashova T.A. , Sukov S.A. Gibridnaya superkomp'yuternaya platforma i razrabotka prilozheniy dlya resheniya zadach mekhaniki sploshnoy sredy setochnymi metodami [Hybrid supercomputer platform and applications programming for the solution of continuous mechanics problems by grid methods] // Izvestiya YuFU . Tekhnicheskie nauki [IZVESTIYA SFedU . ENGINEERING SCIENCES] . 2012 . No. 6 ( 131 ). P. 105 - 115 .

29. Official documentation and manuals on MPI . URL: http://www.mcs.anl.gov/research/projects/mpi/ (date accessed: 14 . 02 . 2016 ).

30. Official documentation and manuals on OpenMP . URL: http://openmp.org/wp/, https://computing.llnl.gov/tutorials/openMP/ (date accessed: 14 . 02 . 2016 ).

31. CUDA Toolkit Documentation v. 7 .5. Santa Clara (CA, USA): NVIDIA Corporation . URL: http://docs.nvidia.com/cuda/index. html#axzz3tx7h9bwP (date accessed: 14.02 . 2016 ).