Численное моделирование акустических волн в жидком кристалле с использованием технологии CUDA Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 22, Специальный выпуск 1, 2017

Численное моделирование акустических волн в жидком кристалле с использованием технологии CUDA

И. В. СмолЕхо*, О. В. Садовская, В. М. Слдовский

Институт вычислительного моделирования СО РАН, Красноярск, Россия *Контактный e-mail: ismol@icm.krasn.ru

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

Ключевые слова: жидкий кристалл, моментная среда, разностная схема, динамика, параллельный вычислительный алгоритм, технология СИБЛ.

Введение

Жидкие кристаллы представляют собой промежуточное агрегатное состояние вещества между твердыми кристаллами и изотропной жидкостью, в котором проявляются одновременно свойства упругости и текучести. Жидкокристаллическая фаза существует в интервале от температуры плавления до некоторой более высокой температуры, при нагреве до которой вещество переходит в обычную жидкость. Ниже этого интервала вещество является твердым кристаллом. Подвижность молекул жидких кристаллов позволяет воздействием внешних сил изменять их ориентацию и таким образом управлять их свойствами. Основные свойства жидких кристаллов описаны в [1]. Впервые модель жидкого кристалла как микрополярной среды с вращательными степенями свободы была предложена Эриксеном [2]. Аналогичные модели изучались в [3, 4] и других работах. В [5, 6] построена упрощенная модель нематического жидкого кристалла как акустической микронеоднородной среды с вращающимися частицами. Параллельный вычислительный алгоритм для решения этой системы представлен в [7, 8].

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

© ИВТ СО РАН, 2017

1. Математическая модель

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

pu,t = -р,х - q,y, pv,t = q,x - р,у, ju,t = 2q + ¡jx,x + ру>у, <p,t = ш, ptt = -k(u,x + viV) + [Tt, q,t = a(y,x - uy) - 2а(ш + q/r¡), (1)

px,t = 1ш>х, py,t = ^ш,у,

cTt = («11 Tx + «12Ty) x + («12Tx + «22 ТД^ - [T{u^ + v,y) +2q2/'q.

Здесь u и v — компоненты вектора скорости; ш — угловая скорость; < — угол поворота молекул кристалла; р — гидростатическое давление; q — касательное напряжение; рх и ру — моментные напряжения; T — абсолютная температура; р — плотность; j — момент инерции; к — модуль объемного сжатия; a — модуль упругого сопротивления вращению; 7 — модуль упругого сопротивления изменению кривизны; r¡ — коэффициент вязкости; [ — коэффициент теплового расширения; с — удельная теплоемкость; «п, «12 и «22 — коэффициенты тензора теплопроводности: «11 = «1 cos2 < + «2 sin2 <, «12 = («1 - «2)sin< cos<, «22 = «1 sin2 < + «2 cos2 <, («1 и «2 — коэффициенты теплопроводности в направлении ориентации молекул жидкого кристалла и в поперечном направлении), нижние индексы t и х, у означают производные по времени и пространственным переменным соответственно.

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

Продифференцируем первое уравнение системы (1) по у, второе — по х. После вычитания получим

p(u,y - víx)j = -Aq,

где A — оператор Лапласа. С учетом уравнений для рх и ру после дифференцирования уравнения вращательного движения получим систему двух уравнений второго порядка для касательного напряжения и угловой скорости ш:

2 a a 2

q,tt +--q,t + 2awit = - Aq, ш,ы -- q,t = - Аш. (2)

Начальные данные для системы (2):

lt=0

q,t\t=o = a(v% - u01 - 2а(ш0 + £) = -2а(ш0 + £)

ш I t=o = ш,

1 2 п0

0 0 0

t \ t=0 ~¡\2<i + рх,Х + ру,у)

где u0, V0, ш0, q0, рРх, р0 — константы, заданные в начальный момент времени.

2. Численный алгоритм

Для решения системы уравнений (2) с начальными данными (1) и граничными условиями разработан вычислительный алгоритм. Искомыми величинами являются касательное напряжение д и угловая скорость ш внутри расчетной области. Граничные условия, моделирующие внешнее воздействие, задаются в терминах р ш, а условия симметрии напряженного состояния жидкого кристалла формулируются в терминах производных и шХ или и ш,у в зависимости от положения линии симметрии. Используется явная конечно-разностная схема "крест" второго порядка точности по х, у и ¿.

Аппроксимируем уравнения системы (2) на каждом временном слое заменой производных по времени и пространству конечными разностями:

п+1 9 п I „п-1 пп+1 пп-1 , ,

$11 ,72 - 2 $11 ,72 + $11 ,72 + а $31,32 - $31,32 + а Ш71 ,72 - Ш7

(△)2 г] Аг Аг

/ Пп _ О Пп I Пп Пп _ О Пп I Пп \

а $31+1,32 2 $31,32 + 431-1,32 + Чл ,32+1 2 $31 ,32 + $31 ,32 - М

Р^ (Ах)2 + (Ау )2 У

. п+1 2 шп + шп-1 1 пп+1 па-1 Ш71 ,72 - 2 Ш71 ,72 + Ш71 ,72 1 ?71 ,72 - ?71 ,72

(А )2 А

7/ ШЛ + 1,Л - 2Ш31,32 + ШЛ-1,^2 + Ш31,32 + 1 - 2 Ш31,32 + Ш31,32-1 \

А (Ах)2 + (Ау )2 )

где ]1 = 2, М1 — 1, ]2 = 2, М2 — 1. Через АЬ и Ах, Ау обозначены шаги равномерной сетки по времени и пространственным переменным; N1 и Ы2 — размерности сетки в направлениях х и у; нижние индексы ]1 и ]2 характеризуют пространственное положение ячейки сетки, верхний индекс (п + 1) соответствует текущему моменту времени, индексы п и (п — 1) — предыдущим моментам.

Далее из второго уравнения выразим :

шП+1 = 2 , ,П Ш™-1 + А I ПП+1 /7га-М +

ш31,32 = 2 Ш31,32 — ш31,32 + — у31,32 — $31,32 ) +

/ л I \2 / ,п О , I , , О , I , \

+ 7 (А) Ш31+1,32 — 2 Ш31,32 + ш31-1,32 + ш31,32 + 1 — 2 ш31,32 + ш 31,32-Л

3 V (Ах)2 (Ау )2 У

После подстановки (4) в первое уравнение системы (2) получим формулу для :

±31,32

(а а 1 А п+1 _ 2 п

Ч + + Щ]2)= Щ)2 +

3 г]Аг (А£)2)(Аг)2 Чзъп а + ___\пп-1 V

2 г]АЬ (Аí)2)+ А* VШJ1,J2 Ш^2)

/ Пп _ о Пп Пп Пп _ о ^п Пп

+ а $31 + 1,32 2 $31,32 + $31-1,32 + $31,32 +1 2 ^1,^2 + ^1,^2 -1

+ Р V (Ах)2 + (Ау)2 ,

Л х / / О , I . ,П . ,П О , I , ,п

а7 А ш31+1,32 — 2ш31,32 + шл-1,^2 , шл,л+1 — 2шл,^2 + шл,^2-

. С 1 ■ 2_„, 1 2_„, 1 2 I Л02 I ^_Ц Ц 2_ЦЦ2 1 |

V (Ах)2 + (А у )2 У

Таким образом, на новом слое по времени нужно сначала решать уравнение (5) для д, а затем уравнение (4) для ш при заданных начальных данных и граничных условиях.

3. Устойчивость схемы

При исследовании устойчивости конечно-разностной схемы для простоты пренебрежем вязким слагаемым, устремив г/ ^ то. Это упрощение основано на том, что вязкость увеличивает запас устойчивости схемы, в чем можно убедиться с помощью численных экспериментов. В соответствии с методом Фурье положим

31,32

Хп q ei(jlal+Í2a2)

ш

31,32

xa q ei(iiai+i2a2)

(6)

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

Л2 - 2 Л +1 4 a /sin2(ai/2) sin2(a2/2) Л2 -

+ y4yAxy + yA)H )q+a^TUJ

(Ai)2 р \ (Дх)2

Л2 - 2 Л +1 4j /sin2(ai/2) sin2(a2/2) \\ (Д)2 + У Ч (Дх)2 + (Ду)2 ))

1 Л2 - 1 ш--т ——— q

0,

0.

3 \ (△ж)2 (Ду )2 )) з △

Составив матрицу из коэффициентов при д и С, запишем следующее характеристическое уравнение

( А - 1)2

4 a

/ Л чо +--ЛА

(Ai)2 р

1 Л2 1

Л2 1

a

Д

(Л - 1)2 + ИлА

(Д )2

0,

. sin2(ai/2) sin2(a2/2) ^

где А = ————--1--————. Введя обозначения

(Д х)2 (Д )2

а = аА (At )2, Ъ=- А (At )2, с= a (At )2, вычислим определитель

(1 + с)(Л2 - 1)2 + 4 Л (Л - 1)2 (а + Ь- 1) + 16 Л2 аЬ = 0.

Далее рассмотрим три случая при разных значениях а, b и с.

1. Если 6 = 0 (т. е. если 7 = 0), то характеристическое уравнение решается в явном виде. Два его корня совпадают и равны 1, еще два являются решениями квадратного уравнения (1 + с)( Л + 1)2 + 4 Л (а - 1) = 0. По теореме Виета произведение этих корней равно 1. Поэтому в случае действительных корней, когда дискриминант положительный, один из них строго больше единицы. Это случай неустойчивости схемы. Схема устойчива, если дискриминант меньше либо равен нулю. В этом случае корни комплексно-сопряженные, следовательно, |Л1| = |Л2| = 1. Таким образом, из условия неположительности дискриминанта

1 - а\2

1 - 2-) - 1 ^ 0

1+

а ^ 1

вытекает спектральное условие устойчивости схемы

a , / 1 1 ^ -1

Р

a (л )2 ^ Í 1 1 \

~о ( t] ^ V(Ax)2 + (Aу)2)

2. Если а = 0 (т. е. если а = 0), то соответствующее квадратное уравнение записывается в виде (1 + с)(Л + 1)2 + 4 Л (Ь — 1) = 0. Повторяя предыдущие рассуждения, получим условие неположительности дискриминанта Ь ^ 1, из которого вытекает условие устойчивости схемы ^

](Л)2 ^ ((дХ^ + (д^) •

3. Если а + Ь = 1, то характеристическое уравнение приводится к биквадратному уравнению (1 + с)( Л2 — 1)2 + 16 Л2 а Ь = 0. Сделав замену г = Л2, получим квадратное уравнение, корни которого комплексно-сопряженные и лежат на единичной окружности при условии неположительности дискриминанта

О 6 \ 2

1 — 8-) — 1 ^ 0 ^ 4а6 ^ 1 + с.

1+

Так как а + Ъ = 1, последнее неравенство приводится к виду 4а6 ^ (а + Ь)2 + с, т.е. 0 ^ (а — 6)2 + с, и автоматически выполняется. Условие а + 6=1 означает, что

а 7\/Л \2 (81п2(а1/2) в1п2(а2/2)\-1

(а+1) )2

р ' V (Лх)2 (Лу )2 / •

Для фиксированных значений а1 и а2 такой выбор шага по времени дает |Л| = 1. Соответствующее решение (6) ограничено. Не вычисляя корней характеристического уравнения, можно доказать, что оно остается ограниченным и при меньших значениях шага по времени.

Докажем более общее утверждение. Пусть ||ига|| — произвольная норма на пространстве решений двухслойной однородной разностной схемы с постоянным оператором шага Л, частным случаем которой является схема (4), (5). Тогда если при переходе на новый слой по времени норма решения не увеличивается в расчетах с шагом ДЬ, то она не увеличивается и в расчетах с меньшим шагом (Д£ < Д£).

Доказательство существенно опирается на свойство выпуклости нормы — следствие ее положительной однородности и неравенства треугольника:

||хи + (1 — х)иЦ ^ х ||и|| + (1 — х)||и||, 0 < х < 1.

Переход на новый слой в схеме осуществляется по формулам

ип+1 =ип + ДгЛип, ип+1 = ип + ДЛи", с помощью которых можно показать, что

ип+1 = ип + х ДЛип = х ип+1 + (1 — х)ип, х = Дt/Д е (0,1). По предположению, ||ига+1|| ^ ||ига||. Следовательно,

|К+1|| ^х ||ига+1|| + (1 — х)|К|| ^ |К||,

что и требовалось доказать.

Из соображения ограниченности решений ип = (с[п, шп) вида (6) для всевозможных значений параметров а1 и а2, получим следующее условие:

'а , 1\,Л^2 / ( 1 + 1 V1

(а+1) )2 < - (7)

Строго говоря, это условие является необходимым спектральным условием устойчивости схемы только в двух частных случаях, рассмотренных в п. 1 и 2. В общем случае оно гарантирует ограниченность решений вида (6), что весьма важно с практической точки зрения, но не известно, являются ли такие решения растущими, если оно нарушено. Кроме того, ни условие (7), ни необходимое спектральное условие устойчивости Фурье, вообще говоря, не являются достаточными и, таким образом, не гарантируют устойчивости, которая может зависеть, например, от способа аппроксимации граничных условий. Однако на практике расчеты по описанной схеме при выборе шага по времени в соответствии с условием (7) показали устойчивые вычислительные результаты.

4. Параллельная программа

Описанный выше алгоритм численного решения системы уравнений (2) по формулам (4), (5) реализован в виде параллельной программы. Программа написана на языке Си с применением технологии CUDA (Compute Unified Device Architecture) [9], позволяющей существенно увеличить вычислительную производительность благодаря использованию графических ускорителей видеокарт.

Вычисления производятся на GPU (Graphical Processing Unit). Расчетная область разбивается на квадратные блоки. Благодаря идентификаторам, имеющимся в CUDA, каждой ячейке конечно-разностной сетки в пределах блока ставится в соответствие нить графического устройства. В параллельном режиме нити выполняют однотипные операции в ячейках по расчету решения на каждом шаге по времени.

Параллельная программа имеет следующую структуру:

• задание размерностей конечно-разностной сетки и всех необходимых констант (на CPU);

• описание одномерных массивов для касательного напряжения и угловой скорости (на CPU);

• задание начальных данных для этих величин в узлах конечно-разностной сетки (на CPU);

• описание событий start и stop для замера времени выполнения программы на GPU, начало замера времени, начало параллельной части работы программы;

• копирование необходимых для расчетов констант из CPU в GPU;

• описание массивов для угловой скорости и касательного напряжения, а также для всех необходимых вспомогательных величин, выделение под них памяти (на GPU);

• копирование данных из CPU в GPU (массивов искомых величин);

• задание переменных типа dim3 для количества блоков в сетке и количества нитей в каждом из таких блоков (на GPU);

• основной расчетный цикл по времени, в котором последовательно выполняются процедуры-ядра (на GPU):

— задание граничных условий;

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

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

Рис. 1. Ускорение GPU по сравнению с CPU

— копирование результатов счета из GPU в CPU (массивов угловой скорости и касательного напряжения) в контрольных точках (на некоторых шагах по времени);

• освобождение памяти переменных на GPU;

• окончание замера времени работы на GPU, печать этого времени, уничтожение событий start и stop, завершение параллельной части работы программы;

• освобождение памяти переменных на CPU, завершение работы программы. Эффективность параллельного кода для решения исходной системы уравнений (1),

описывающей динамику жидкого кристалла, исследовалась в [7]. Было проведено большое количество расчетов при разных размерностях сетки. Сравнивалось время счета двух версий программы — параллельной и последовательной. На рис. 1 показана зависимость ускорения параллельной программы от размерности сетки N х N конечно-разностной схемы. Здесь N принимает значения: 10, 100, 200, 400, 600, ..., 3000, 3200. Получено ускорение работы параллельной программы примерно в 25 раз по сравнению с последовательной версией при размерностях сетки 1000 х 1000 и выше.

5. Точное решение

Для одномерной задачи о действии касательного напряжения на границе проводилось сравнение численного решения с помощью описанной параллельной программы и точного решения. Подставим q = £ ег(?1-кУ), . = Ш ег(^г-ку") в уравнения системы (2). После упрощения получим

п2 2га / ак2\ 2г/* ^ /7к2 п2<

2 2 а а к 2 к 2 -/2 +-- +-)д + 2га/Ш = 0,--- д + ( — — /2)Ш = 0.

Вычислив определитель этой системы, найдем выражение для к±:

к ± =

\

Щл ±< Id? - 4 ^ (f 2 - ^ - ^ )\ d = + l)f - 2г

у pj V ц j J) \р jj nj

Здесь / — частота, к± = к± + г к± — волновые числа.

На рис. 2 и 3 представлены характерные дисперсионные кривые: зависимость фазовой скорости с± = //Яе к± от частоты у = //(2^) в герцах и зависимость декремента затухания А± = — 1/1тк± от частоты у (для к + и к- — слева и справа соответственно).

[м/с] а

0 200 400 600 800 1000 и [Д

с [км/с]

А+ [мкм]

Рис. 2. Зависимость фазовой скорости от частоты: а — для к+, б — для к

а б

Л~ [мкм]

400 350 300 250 200 150

0 200 400 600 800 1000 и [Ь

25 20 15 10

5

0 200 400 600 800 1000 и [Г

Рис. 3. Зависимость декремента затухания от частоты: а — для к+, б — для к

ш [ГПа]

1 0.8 0.6 0.4 0.2 0

0 0.5 1 1.5

2.5 3 3.5 у [мкм] 0 0.5 1 1.5 2 2.5 3 3.5 у [мкм] Рис. 4. Зависимость угловой скорости от координаты: а — для к+, б — для к-

а

Пунктирная линия соответствует собственной частоте вращательного движения частиц кристалла: и* = [р.Ц) ¡ж.

Вычисления проведены для жидкого кристалла 5ЦБ со следующими параметрами: р = 1022 кг/м3, ^ = 1.33 ■ 10-10 кг/м, а = 0.161 ГПа, 7 = 10 мкН, ^ =10 Па ■ с. Для этого кристалла и* = 350 МГц. Размер расчетной области 4 мкм.

Начальные данные и граничные условия для одномерной задачи определяются в соответствии с точным решением: Re q = ек2У (giCOs( f t — к\у) — q2 sin( f t — кiy)), Re ш = ек2У (wi cos(/t — k\у) — ш2 sin(/t — кiy)).

На рис. 4 показаны результаты численного решения с описанными выше параметрами: зависимость ш от у для к+ и к- в один из моментов времени. Размерность сетки составляет 1000 ячеек. Относительная погрешность равна 3 • 10-3 в расчетах для к+, 5 • 10-4 — для к-.

6. Результаты расчетов

Для демонстрации работы программы проведена серия расчетов на высокопроизводительном вычислительном сервере Flagman с графическими вычислителями Tesla C2050 Института вычислительного моделирования СО РАН.

На рис. 5 представлены результаты расчетов для задачи о действии сосредоточенного касательного напряжения q = q8(x)8(t) на верхней границе расчетной области. На этой границе q = д, если |x — xc| ^ Ax и t ^ At, и q = 0 в противном случае, ш,у = 0. Нижняя граница закреплена: q = 0, ш = 0. На правой и левой границах заданы условия периодичности. В начальный момент времени касательное напряжение q и угловая скорость ш равны нулю. Вычисления проведены на прямоугольной области 10 х 4 мкм. Нагрузка приложена в точке xc = 5 мкм. Расчеты выполнены для жидкого кристалла 5ЦБ, параметры которого приведены в предыдущем параграфе. Размерность конечно-разностной сетки 2560 х 1024 ячеек.

На рис. 6 показаны результаты расчетов для задачи о действии шести П-образных импульсов касательного напряжения на боковых границах расчетной области (по три с каждой стороны). Время действия каждого импульса и промежуток между ними 8 нс. Граничные условия на левой и правой границах: q = д, если |у — yd ^ I, и q = 0 в противном случае, ш,х = 0. На верхней и нижней границах заданы условия периодичности. Вычисления проведены на области 100 х 40 мкм. Центр зоны, в которой действует нагрузка, ус = 20 мкм, радиус этой зоны I = 10 мкм. Размерность сетки 2560 х 1024 ячеек.

a

Рис. 5. Действие сосредоточенного касательного напряжения на верхней границе: линии уровня касательного напряжения (а) и угловой скорости (б) (1000-й и 2000-й шаги по времени)

Рис. 6. Действие П-образных импульсов касательного напряжения на боковых границах: линии уровня касательного напряжения (1000-й, 2000-й, 3000-й и 3500-й шаги по времени)

Рис. 7. Действие постоянного П-образного импульса касательного напряжения на верхней границе: поверхности уровня угловой скорости (1000-й и 2000-й шаги по времени)

Рис. 8. Периодическое действие касательного напряжения на верхней границе: поверхности уровня угловой скорости (1000-й и 2000-й шаги по времени)

В отличие от предыдущей задачи, здесь взяты параметры j = 1.33 • 10 7 кг/м, 7 =1 мН, •п = 100 Па • с.

Поверхности уровня угловой скорости для задачи о действии постоянного П-образно-го импульса касательного напряжения на части границы изображены на рис. 7. На верхней границе q = q при \х — хс\ ^ /, ш,у = 0. На правой и левой границах заданы условия периодичности, нижняя граница закреплена. Расчетная область 100 х 40 мкм, хс = 50 мкм, I = 25 мкм.

На рис. 8 представлены результаты расчетов для задачи о периодическом действии касательного напряжения на части верхней границы. На этой границе q = q sin(2nut), если \х — хс\ ^ /,и q = 0 в противном случае, ш,у = 0. Нижняя граница закреплена, на левой и правой границах заданы условия периодичности. В расчетах хс = 5 мкм, I = 2.5 мкм, размеры области 10 х 4мкм. Частота v равна собственной частоте v* = 350 МГц вращательного движения кристалла 5ЦБ. Результаты расчетов показывают рост амплитуды угловой скорости со временем, характерный для резонансного возбуждения среды.

Заключение

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

Благодарности. Работа выполнена при финансовой поддержке Комплексной программы фундаментальных исследований СО РАН № 11.2П "Интеграция и развитие" (проект № 0356-2016-0728).

Список литературы / References

[1] Blinov, L.M. Structure and properties of liquid crystals. Heidelberg; New York; Dordrecht; London: Springer, 2011. 440 p.

[2] Ericksen, J.L. Conservation laws for liquid crystals // Trans. Soc. Rheol. 1961. Vol. 5. P. 23-34.

[3] Leslie, F.M. Some constitutive equations for liquid crystals // Arch. Ration. Mech. Anal. 1968. Vol. 28. P. 265-283.

[4] Аэро Э.Л., Булыгин А.Н. Уравнения движения нематических жидких кристаллов // Прикл. математика и механика. 1971. Т. 35, вып. 5. С. 879-891.

Aero, E.L., Bulygin, A.N. Equations of motion of nematic liquid-crystal media //J. Appl. Math. Mech. 1971. Vol. 35, No. 5. P. 831-843.

[5] Садовский В.М., Садовская О.В. Об акустическом приближении термомеханической модели жидкого кристалла // Физическая мезомаханика. 2013. Т. 16, № 3. С. 55-62. Sadovskii, V.M., Sadovskaya, O.V. On the acoustic approximation of thermomechanical description of a liquid crystal // Phys. Mesomech. 2013. Vol. 16, No. 4. P. 312-318.

[6] Sadovskii, V.M. Equations of the dynamics of a liquid crystal under the influence of weak mechanical and thermal perturbations // AIP Conf. Proc. 2014. Vol. 1629. P. 311-318.

[7]

[8]

[9]

Поступила в 'редакцию 18 января 2017 г.

Numerical modeling of acoustic waves in a liquid crystal using CUDA technology

Smolekho, Irina V.*, Sadovskaya, Oxana V., Sadovskii, Vladimir M.

Institute of Computational Modeling SB RAS, Krasnoyarsk, 660036, Russia * Corresponding author: Smolekho, Irina V., e-mail: ismol@icm.krasn.ru

One of the approaches to mathematical modeling for deformation of a liquid crystal is based on the assumption that a liquid crystal is a fine-dispersed continuum, at each point of which the elongated particles — molecules or domains of co-oriented molecules — can move according to laws of the viscous fluid dynamics and can rotate relative to a fluid. In this approach the Cosserat continuum theory is applicable, where along with translational motion, characterized by the velocity vector, the rotational degrees of freedom for particles with the angular velocity vector are considered and, along with the stress tensor with nonsymmetrical components, the nonsymmetrical tensor of couple stresses is introduced. Within the framework of acoustic approximation of the model, which describes thermomechanical behavior of a liquid crystal taking into account the couple stresses, the system of two differential equations of the second-order was obtained for the tangential stress and angular velocity in two-dimensional case.

The algorithm for numerical solution of the obtained system of equations under given initial data and boundary conditions is considered. This system is solved using the explicit finite-difference scheme "cross" of the second-order approximation. The stability condition for the scheme is obtained. The algorithm is realized as a parallel program in the C language using the CUDA technology for computing systems with graphics accelerators. A series of computations of acoustic waves in a liquid crystal was performed to demonstrate the efficiency of proposed program. The comparison of numerical solution and the exact solution was carried out.

Keywords: liquid crystal, couple-stress medium, dynamics, difference scheme, parallel computational algorithm, CUDA technology.

Acknowledgements. This research was supported by the Complex Fundamental Research Program No. II.2P "Integration and Development" of SB RAS (project No. 03562016-0728).

Received 18 January 2017

Sadovskaya, O.V. Numerical simulation of the dynamics of a liquid crystal in the case of

plane strain using GPUs // AIP Conf. Proc. 2014. Vol. 1629. P. 303-310.

Смолехо И.В. Параллельная реализация алгоритма для описания термоупругих волн в

жидких кристаллах // Молодой ученый. 2015. № 11(91). С. 107-112.

Smolekho, I.V. Parallel implementation of the algorithm for description of thermoelastic

waves in liquid crystals // Molodoy Uchenyy. 2015. No. 11(91). P. 107-112. (In Russ.)

Farber, R. CUDA application design and development. Amsterdam; Boston; Tokyo: Morgan Kaufmann / Elsevier, 2011. 315 p.

© ICT SB RAS, 2017