=Paper=
{{Paper
|id=Vol-1839/MIT2016-p40
|storemode=property
|title= Numerical analysis of acoustic waves in a liquid crystal taking into account couple-stress interaction
|pdfUrl=https://ceur-ws.org/Vol-1839/MIT2016-p40.pdf
|volume=Vol-1839
|authors=Irina Smolekho,Oxana Sadovskaya,Vladimir Sadovskii
}}
== Numerical analysis of acoustic waves in a liquid crystal taking into account couple-stress interaction==
https://ceur-ws.org/Vol-1839/MIT2016-p40.pdf
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Numerical Analysis of Acoustic Waves
in a Liquid Crystal Taking Into Account
Couple-Stress Interaction
Irina Smolekho, Oxana Sadovskaya, and Vladimir Sadovskii
Institute of Computational Modeling SB RAS,
Akademgorodok 50/44, 660036 Krasnoyarsk, Russia
{ismol,o_sadov,sadov}@icm.krasn.ru
http://icm.krasn.ru
Abstract. Based on the mathematical model describing the thermo-
mechanical behavior of a liquid crystal, which takes into account the
couple-stress interactions, the system of two differential equations for
tangential stress and angular velocity was obtained in two-dimensional
case. Computational algorithm for numerical solution of this system of
equations of the second order under given initial data and boundary
conditions is worked out. The algorithm is implemented as a parallel
program in the C language using the CUDA technology for computer
systems with graphic accelerators. A series of numerical computations of
acoustical waves in a liquid crystal was carried out to demonstrate the
efficiency of proposed parallel program.
Keywords: liquid crystal, couple-stress medium, dynamics, finite-diffe-
rence scheme, parallel computational algorithm, CUDA technology
1 Introduction
Liquid crystal – it is an intermediate state of matter, which appears at the same
time the properties of elasticity and fluidity. The liquid crystal phase is formed
during melting of a number of organic substances. It exists in a range from
the melting temperature to a higher temperature, when heated to a substance
which becomes an ordinary liquid. Below this range the substance is a solid
crystal. In the liquid crystal state may be some organic compounds consisting
of molecules of an elongated shape (in the form of elongated rods or plates)
having parallel folding of such molecules. Liquid crystal as a liquid can take
the container shape in which it is placed. However, apart from this property,
combining it with a liquid, it has the property, characteristic of crystals, – the
presence of the order of the spatial orientation of the molecules. The liquid
crystal molecules are oriented in a direction which is determined by the unit
vector, called “director”. Depending on the type of ordering axes of liquid crystal
molecules, they are divided into three types: nematic, smectic and cholesteric.
Nematic and smectic liquid crystals are characterized by a parallel arrangement
473
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
of the molecules. Cholesteric liquid crystals are a kind of nematic liquid crystals,
but they lack long-range order. In this paper we consider a nematic liquid crystal.
Nematic liquid crystals are characterized by the orientation of the longitudinal
axes of the molecules along a certain direction (long-range orientational order
is characteristic for them). Molecules continuously slide in the direction of their
long axes, revolving around them, but at the same time retain orientational
order: the long axes are directed along a preferred direction. In a nematic state,
not all molecules have the same orientation. Because of the lack of fixing of
molecules the orientation of the director changes. Since the director at different
sections is oriented differently, the regions with different directions of director –
the domains – appear in a liquid crystal. At the interfaces of domains the light
refractive index changes, so liquid crystals become hazy. The main properties of
liquid crystals are described in [1].
Liquid crystals find many applications because of strong dependence of their
properties on external influences. Due to their ability to reflect light, the liquid
crystals are widely used in laboratories and in the art as convenient means of
visualizing the thermal fields and temperature changes. Due to the high reso-
lution of the liquid crystals, they can be used in microelectronics technology to
detect defects in chips, that improves reliability of the circuit. With the help of
liquid crystals in medicine can directly observe the distribution of the human
body surface temperature, and it is important to identify foci of inflammation
hidden under the skin. Liquid crystals found the most widely usage in digital
technology. Currently, color LCD screens have even greater range of applications:
mobile phones, personal computers and televisions, which have a small thickness,
low power, high resolution and brightness.
One of the approaches to the construction of a mathematical model to de-
scribe the behavior of liquid crystals is based on the representation of a liquid
crystal medium as a fine-dispersed continuum. At each point of this continuum,
the domains of a liquid crystal can move in accordance with laws of the dyna-
mics of viscous or inviscid liquid and can rotate relative to a liquid, encounte-
ring resistance to rotation. The models of liquid crystals have been proposed by
Eriksen [2], Leslie [3], Aero [4] and other authors. This paper is devoted to the
numerical solution of differential equations of the second order for tangential
stress and angular velocity, obtained from the system of equations describing
the thermomechanical behavior of a liquid crystal in the two-dimensional case.
2 Governing Equations
In the framework of acoustic approximation, the mathematical model of a liquid
crystal without taking into account the couple stresses is described in [5, 6].
The system of equations of this model includes the equations of translational
and rotational motion, the equation for the angle of rotation, the constitutive
equations for pressure and tangential stress, as well as the equation of anisotropic
heat conduction with variable coefficients. Parallel computational algorithm for
the solution of this system is represented in [7, 8].
474
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
In two-dimensional case the complete system of equations describing the be-
havior of a liquid crystal under weak acoustic perturbations taking into account
the couple stresses is as follows:
ρ u,t = −p,x − q,y , ρ v,t = q,x − p,y ,
j ω,t = 2 q + µx,x + µy,y , φ,t = ω,
( ) ( ) ( )
p,t = −k u,x + v,y + β T,t , q,t = α v,x − u,y − 2 α ω + q/η ,
(1)
µx,t = γ ω,x , µy,t = γ ω,y ,
( ) ( )
c T,t = æ11 T,x + æ12 T,y ,x + æ12 T,x + æ22 T,y ,y −
( )
− β T u,x + v,y + 2 q 2 /η.
Here u and v are the projections of the velocity vector on the coordinate axes,
ω and φ are the angular velocity and the rotation angle, p is the hydrostatic
pressure, q is the tangential stress, µx and µy are the couple stresses, T is
the absolute temperature, ρ is the density, j is the moment of inertia, k is
the bulk compression modulus, α is the modulus of elastic resistance to rota-
tion, η is the viscosity coefficient, c and β are the coefficients of heat capacity
and thermal expansion, æ11 , æ12 and æ22 are the components of the thermal
conductivity tensor: æ11 = æ1 cos2 φ + æ2 sin2 φ, æ12 = (æ1 − æ2 ) sin φ cos φ,
æ22 = æ1 sin2 φ + æ2 cos2 φ, (æ1 and æ2 are the thermal conductivity coef-
ficients of a liquid crystal in the direction of molecular orientation and in the
transverse direction). Subscripts after a comma denote the partial derivatives
with respect to time t and spatial variables x and y.
The system (1) includes the equations of translational and rotational motion,
the equation for the angle of rotation, the equations of state for pressure and tan-
gential stress, the equations for couple stresses, and the equation of anisotropic
heat conduction with variable coefficients.
Let’s consider how to obtain the system of equations of the second order for
tangential stress and angular velocity. Differentiating the first equation of (1)
by x, the second equation by y and subtracting the second from the first, we
find:
( )
ρ u,y − v,x ,t = −△q,
where △ is the Laplace operator. In view of this expression and also expressions
for µx,t and µy,t , after differentiation of corresponding equations of the system (1)
by t, we obtain a separate subsystem for the tangential stress q and the angular
velocity ω:
2α α
q,tt + q,t + 2 α ω,t = △q,
η ρ
(2)
2 γ
ω,tt − q,t = △ω.
j j
475
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Initial data for the system (2) have the following form:
( 0 ) ( q0 ) ( q0 )
q t=0 = q 0 , q,t t=0 = α v,x − u0,y − 2 α ω 0 + = −2 α ω 0 + ,
η η
(3)
1( ) 2 q0
ω t=0 = ω 0 , ω,t t=0 = 2 q 0 + µ0x,x + µ0y,y = ,
j j
where u0 , v 0 , ω 0 , q 0 , µ0x , µ0y are given constants at the initial time moment.
The fourth-order equation for q can be derived from the subsystem (2). So,
expressing ω,t from the first equation of (2) and substituting it into the second
equation, differentiated by t, we obtain the next chain of equations:
1 (α 2α ) 2 γ
ω,t = △q − q,tt − q,t , ω,ttt = q,tt + △ω,t ,
2α ρ η j j
2α 4α (α γ ) 2αγ α γ 2
q,tttt + q,ttt + q,tt − + △q,tt − △q,t = − △ q.
η j ρ j jη ρ j
Initial data for the corresponding Cauchy problem are as follows:
( 0 ) ( q0 ) ( q0 )
q t=0 = q 0 , q,t t=0 = α v,x − u0,y − 2 α ω 0 + = −2 α ω 0 + ,
η η
α 2 α( 0 ) 2α 0 [α (α 1 ) 0]
q,tt t=0 = △q 0 − 2 q + µ0x,x + µ0y,y − q,t = 4 α ω0 + 2 − q ,
ρ j η η η j
α 0 2α ( 0 ) 2α 0
q,ttt t=0 = △q,t − 2 q,t + γ △ω 0 − q =
ρ j η ,tt
[( 1 α) 1(1 1 α) ]
= 8 α2 − 2 ω0 + + − 2 q0 .
j η η j jη η
3 Finite-Difference Scheme
Computational algorithm is developed for numerical solution of the system of two
second-order equations (2) with the initial data (3). The unknown variables are
the tangential stress q and the angular velocity ω within computational domain.
Boundary conditions are defined in terms of q, ω and also q,x , ω,x , q,y , ω,y .
The explicit finite-difference scheme “cross” of the second order approximation
by x, y and t is used [9].
Equations of the system (2) at each time step are approximated by replacing
the derivatives with respect to time and spatial variables by the finite differences:
qjn+1
1 ,j2
− 2 qjn1 ,j2 + qjn−1
1 ,j2
n+1 n−1
α qj1 ,j2 − qj1 ,j2 ωjn+1
1 ,j2
− ωjn−1
1 ,j2
+ + α =
(△t)2 η △t △t
( )
α qjn1 +1,j2 − 2 qjn1 ,j2 + qjn1 −1,j2 qjn1 ,j2 +1 − 2 qjn1 ,j2 + qjn1 ,j2 −1
= + ,
ρ (△x)2 (△y)2
476
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
ωjn+1
1 ,j2
− 2 ωjn1 ,j2 + ωjn−1
1 ,j2
n+1 n−1
1 qj1 ,j2 − qj1 ,j2
− =
(△t)2 j △t
( )
γ ωjn1 +1,j2 − 2 ωjn1 ,j2 + ωjn1 −1,j2 ωjn1 ,j2 +1 − 2 ωjn1 ,j2 + ωjn1 ,j2 −1
= + ,
j (△x)2 (△y)2
where j1 = 2, N1 − 1, j2 = 2, N2 − 1. Next, ωjn+1 1 ,j2
can be expressed from the
second equation:
( )
△t n+1
ωjn+1
1 ,j2
= 2 ω n
j1 ,j2 − ω n−1
j1 ,j2 + qj1 ,j2 − q n−1
j1 ,j2 +
j
2( ) (4)
γ (△t) ωjn1 +1,j2 − 2 ωjn1 ,j2 + ωjn1 −1,j2 ωjn1 ,j2 +1 − 2 ωjn1 ,j2 + ωjn1 ,j2 −1
+ 2 + 2 .
j (△x) (△y)
Substituting (4) into the first equation, we obtain the formula for qjn+1
1 ,j2
:
(α α 1 ) n+1 2
+ + q = qn +
j η △t (△t)2 j1 ,j2 (△t)2 j1 ,j2
(α ( )
α 1 ) n−1 2α n−1 n
+ + − q + ωj1 ,j2 − ωj1 ,j2 +
j η △t (△t)2 j1 ,j2 △t
( n ) (5)
α qj1 +1,j2 − 2 qjn1 ,j2 + qjn1 −1,j2 qjn1 ,j2 +1 − 2 qjn1 ,j2 + qjn1 ,j2 −1
+ + −
ρ (△x)2 (△y)2
( n )
α γ △t ωj1 +1,j2 − 2 ωjn1 ,j2 + ωjn1 −1,j2 ωjn1 ,j2 +1 − 2 ωjn1 ,j2 + ωjn1 ,j2 −1
− + .
j (△x)2 (△y)2
Calculating tangential stress by the formula (5) and then angular velocity by the
formula (4) at each time step, one can find numerical solution of the problem.
Finite-difference scheme has the second-order approximation by time and
spatial variables. According to the Lax theorem, the sequence of approximate
solutions converges to the exact solution with the second order, too.
4 Stability of the Scheme
Under analysis of the stability of the finite-difference scheme for simplicity let’s
neglect the viscous term tending η → ∞. This simplification is based on the
assumption that viscosity increases the reserve of stability of the scheme. Ac-
cording to the Fourier method, let
qjn1 ,j2 = λn q̂ ei(j1 α1 +j2 α2 ) , ωjn1 ,j2 = λn ω̂ ei(j1 α1 +j2 α2 ) .
Substituting these values into the first equation of the system (2) and dividing
both sides of the equation by λn ei(j1 α1 +j2 α2 ) , we get:
( )
λ − 2 + 1/λ λ − 1/λ α eiα1 − 2 + e−iα1 eiα2 − 2 + e−iα2
q̂ + α ω̂ = + q̂.
(△t)2 △t ρ (△x)2 (△y)2
477
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Consequently,
( )
λ2 − 2 λ + 1 4 α ( sin2 (α1 /2) sin2 (α2 /2) ) λ2 − 1
+ λ + q̂ + α ω̂ = 0.
(△t)2 ρ (△x)2 (△y)2 △t
After similar calculations for the second equation of (2) we find:
( )
λ2 − 2 λ + 1 4 γ ( sin2 (α1 /2) sin2 (α2 /2) ) 1 λ2 − 1
2
+ λ 2
+ 2
ω̂ − q̂ = 0.
(△t) j (△x) (△y) j △t
To obtain the characteristic equation, we form the matrix of coefficients under
q̂ and ω̂:
(λ − 1)2 4α (λ − 1)(λ + 1)
+ λA α
(△t)2 ρ △t
2 2 2
= 0,
1 λ −1 (λ − 1) 4γ
− + λA
j △t (△t)2 j
sin2 (α1 /2) sin2 (α2 /2)
where A = + . Introducing the notations
(△x)2 (△y)2
α γ α
a= A (△t)2 , b= A (△t)2 , c= (△t)2 ,
ρ j j
one can calculate the determinant:
(λ − 1)4 + 4 λ (λ − 1)2 (a + b) + 16 λ2 a b + (λ2 − 1)2 c = 0,
(1 + c)(λ2 − 1)2 + 4 λ (λ − 1)2 (a + b − 1) + 16 λ2 a b = 0.
So, let us consider three cases with different values a, b and c:
( 1 − a)
1) If b = 0, then (1+c)(λ+1)2 +4 λ (a−1) = 0, λ2 +2 λ 1−2 +1 = 0,
1+c
( 1 − a )2
λ1 = λ̄2 , |λ1 | = |λ2 | = 1, 1−2 − 1 6 0,
1+c
1−a 1−a
−1 6 1 − 2 6 1, 2>2 > 0, a 6 1.
1+c 1+c
Therefore, in this case we obtain the stability condition ∀ α1 , α2 :
α ( 1 1 )−1
(△t)2 6 + .
ρ (△x)2 (△y)2
2) If a = 0, then (1 + c)(λ + 1)2 + 4 λ (b − 1) = 0 and, similarly to the previous
case, one can find that b 6 1. The stability condition is as follows:
γ ( 1 1 )−1
(△t)2 6 2
+ .
j (△x) (△y)2
478
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
3) If a + b = 1, then (1 + c)(λ2 − 1)2 + 16 λ2 a b = 0. Making the substitution
( ab )
z = λ2 , we obtain z 2 − 2 z 1 − 8 + 1 = 0,
1+c
( a b )2 ab
|z1 | = |z2 | = 1, 1−8 − 1 6 0, −1 6 1 − 8 6 1,
1+c 1+c
4 a b 6 1 + c, 4 a b 6 (a + b)2 + c, 0 6 (a − b)2 + c.
The latter condition is satisfied automatically. The condition a + b = 1 means
that
(α γ ) ( sin2 (α /2) sin2 (α /2) )−1
1 2
+ (△t)2 6 + .
ρ j (△x)2 (△y)2
Under such choice of time step, |λ| = 1 for given values α1 and α2 .
Thus, in the general case, the following stability condition takes place:
(α γ) ( 1 1 )−1
+ (△t)2 6 + .
ρ j (△x)2 (△y)2
5 Parallel Program
The described algorithm for numerical solution of the system (2) by the for-
mulas (4), (5) is implemented as a parallel program in the C language using
the CUDA technology for computer systems with graphic accelerators, which
allows to significantly increase the computing performance [10]. GPU (Graphics
Processing Unit) is focused on the implementation of programs with a large
amount of computation. Due to the large number of parallel working cores, it
turns an ordinary computer into a supercomputer with the computing speed of
hundreds of times higher than the PC, using only the computing power of the
CPU. All computations are performed on the GPU, which is a coprocessor to
the CPU. The computational domain is divided into square blocks containing
the same number of threads. Each block is an independent set of interacting
threads, threads of different blocks can not communicate with each other. Due
to the identifiers available in the CUDA, each thread is associated with the mesh
of finite-difference grid. In parallel mode, the threads of a graphic device perform
operations of the same type in the meshes of grid on the calculation of solution
at each time step.
Parallel program has the following structure:
1. Setting the dimensions of finite-difference grid and all the constants used
(on the CPU).
2. Description of one-dimensional arrays for tangential stress and angular ve-
locity (on the CPU).
3. Setting the initial data for these variables at the nodes of the finite-difference
grid (on the CPU).
479
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
4. Description of the events start and stop measuring the program execution
time on the GPU, beginning of the measuring time, beginning of the parallel
part of the program.
5. Copy of the constants needed for computations from the CPU to the GPU.
6. Description of the arrays for angular velocity and tangential stress, and also
for all necessary auxiliary quantities, allocation of memory for them (on the
GPU).
7. Copy of the data from the CPU to the GPU (arrays of the unknown quan-
tities).
8. Setting the variables of the dim3 type for the number of blocks in the grid
and the number of threads in each of these blocks (on GPU).
9. The main computational cycle with respect to time, in which the procedures-
kernels are executed sequentially (on GPU):
(a) setting the boundary conditions in the x direction;
(b) setting the boundary conditions in the y direction;
(c) solving the system of equations for tangential stress and angular velocity
by means of the finite-difference scheme “cross”;
after performing of cores, the barrier synchronization is necessary, to en-
sure completion of the computations by each thread before starting the
next computations;
(d) copy of computational results from the GPU to the CPU (arrays of an-
gular velocity and tangential stress) at the control points (in certain time
steps).
10. Free of memory of the variables on the GPU.
11. Ending of the measuring time on the GPU, print of this time, destruction of
the events start and stop, completion of the parallel part of the program.
12. Free of memory of variables on the CPU, completion of work of the program.
Here you can see the part of the program code for computation of tangential
stress and angular velocity at the internal nodes of finite-difference grid by each
thread of the GPU:
__global__ void syst_qw(int it, double *q, double *q1, \
double *q2, double *w, double *w1, double *w2)
{
int ix,iy,id,idm_x1,idp_x1,idm_x2,idp_x2;
double c1,c2,c3,c4,c5,c6,c7,c8,Delta_q,Delta_w;
ix=threadIdx.x+blockIdx.x*blockDim.x;
iy=threadIdx.y+blockIdx.y*blockDim.y;
if ((ix > 0) && (ix < Nx1Dev-1) && (iy > 0) && (iy < Nx2Dev-1))
{
id=IDX1X2(ix,iy,Nx1Dev);
480
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
idm_x1=IDX1X2(ix-1,iy,Nx1Dev); idp_x1=IDX1X2(ix+1,iy,Nx1Dev);
idm_x2=IDX1X2(ix,iy-1,Nx1Dev); idp_x2=IDX1X2(ix,iy+1,Nx1Dev);
c1=alfaDev/jiDev; c2=gamDev/jiDev;
c3=alfaDev/roDev; c4=1./tauDev/tauDev;
c5=2.*alfaDev/tauDev; c6=-c1*gamDev*tauDev;
c7=tauDev/jiDev; c8=alfaDev/etaDev/tauDev;
Delta_q=(q1[idp_x1]-2.0*q1[id]+q1[idm_x1])/h1Dev/h1Dev;
Delta_q+=(q1[idp_x2]-2.0*q1[id]+q1[idm_x2])/h2Dev/h2Dev;
Delta_w=(w1[idp_x1]-2.0*w1[id]+w1[idm_x1])/h1Dev/h1Dev;
Delta_w+=(w1[idp_x2]-2.0*w1[id]+w1[idm_x2])/h2Dev/h2Dev;
q2[id]=((c1-c4+c8)*q[id]+2.*c4*q1[id]+c5*(w[id]-w1[id])
+c6*Delta_w+c3*Delta_q)/(c1+c4+c8);
w2[id]=2.*w1[id]-w[id]+c7*(q2[id]-q[id])+c2*Delta_w/c4;
__syncthreads();
}
}
__host__ void system_qw(int it, dim3 blocks, dim3 threads, \
dim3 blocks1, dim3 threads1, dim3 blocks2, \
dim3 threads2, double *t, double *q, double *q1, \
double *q2, double *w, double *w1, double *w2)
{
...
syst_qw <<>>(it,q,q1,q2,w,w1,w2);
cudaThreadSynchronize();
}
Previously, in solving the problem without taking into account the couple-
stress interactions, the efficiency of the parallel program for complete system of
equations of a model was analyzed [7]. To evaluate the efficiency of paralleliza-
tion, a large number of computations was performed at different grid dimensions.
The computation time for parallel program and sequential program was com-
pared. Fig. 1 shows a graph of the dependence of acceleration of the parallel
30
25
tCP U /tGP U
20
15
10
5
0 400 800 1200 1600 2000 2400 2800 3200
N
Fig. 1. Acceleration of the program on GPU as compared with CPU
481
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
program on the dimension N × N of a finite-difference grid. Here N takes the
values: 10, 100, 200, 400, 600, ..., 3000, 3200. The acceleration tCP U /tGP U of
the parallel program as compared with the corresponding sequential program is
about 25 times on the grids of dimension 1000 × 1000 and above.
6 Exact Solution
For one-dimensional problem on the action of tangential stress q = q̂ ei(f t−ky) at
one of the boundaries of computational domain, the comparison of the numerical
solution by described parallel program and the exact solution was carried out.
Substituting q = q̂ ei(f t−ky) , ω = ω̂ ei(f t−ky) in the equations of system (2)
after simplification we obtain:
( 2iαf α k2 ) 2if ( γ k2 )
−f 2 + + q̂ + 2 i α f ω̂ = 0, − q̂ + − f 2 ω̂ = 0.
η ρ j j
Calculating the determinant of this system, one can find the expression for k ± :
v √
u
uρj f ( αγ ( 2 2iαf 4 α)
) (α γ ) αγ
±
k = t 2
d± d −4 f − − , d= + f −2 i .
2αγ ρj η j ρ j ηj
Here f is the frequency, k ± = k1± + i k2± are the wave numbers.
The characteristic dispersion curves are represented in Figs. 2 and 3: de-
f f
pendence of the phase velocity c± = ±
on the frequency ν = and
Re k 2π
1
dependence of the damping decrement λ± = − on the frequency ν (for
Im k ±
+ −
k and k – left and right, respectively). The dashed line corresponds√ to the
1 α
eigenfrequency of rotational motion of the particles of a crystal: ν∗ = .
π j
Computations were performed for the liquid crystal 5CB with the next pa-
rameters [11]: ρ = 1022 kg/m3 , j = 1.33·10−10 kg/m, α = 0.161 GPa, γ = 10 µN,
η = 10 Pa · s. For this crystal ν∗ = 350 MHz. The size of a domain is 4 µm.
Initial data and boundary conditions for one-dimensional problem are deter-
mined from the next equations:
( )
Re q̂ = ek2 y q̂1 cos(f t − k1 y) − q̂2 sin(f t − k1 y) ,
( )
Re ω̂ = ek2 y ω̂1 cos(f t − k1 y) − ω̂2 sin(f t − k1 y) .
Figure 4 shows results of numerical solution for described above parameters:
dependence of Re ω on y for k + and k − at one of the instants of time. The
dimension of a finite-difference grid is 1000 meshes. The relative error is 3 · 10−3
in calculations for k + and 5 · 10−4 in calculations for k − .
482
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
a b
300 60
250 50
c− , km/s
200 40
c+ , m/s
150 30
100 20
50 10
0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200
ν, MHz ν, MHz
Fig. 2. Dependence of phase velocity on frequency: a) for k+ , b) for k− .
a b
450 30
400 25
350
20
λ− , µm
λ+ , µm
300
15
250
10
200
150 5
0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200
ν, MHz ν, MHz
Fig. 3. Dependence of damping decrement on frequency: a) for k+ , b) for k− .
a b
2 1
0.8
1
ω, GPa
ω, GPa
0.6
0
0.4
−1
0.2
−2 0
0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4
y, µm y, µm
Fig. 4. Dependence of angular velocity on coordinate: a) for k+ , b) for k− .
483
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
7 Results of Computations
A series of numerical calculations was carried out on the high-performance com-
putational server Flagman with eight graphic solvers Tesla C2050 (448 CUDA
cores on each GPU) of the Institute of Computational Modeling SB RAS to
demonstrate the efficiency of proposed parallel program.
In Fig. 5 one can see the results of computations for the problem on the
action of three Π-shaped impulses of tangential stress on the parts of lateral
boundaries of computational domain. Duration of each impulse and the interval
between them are 8 ns. Initial data are zero. The boundary conditions at the left
and right boundaries: q = q̄, if |y − yc | 6 l, and q = 0, if |y − yc | > l; ω,x = 0.
Here yc is the center of zone, where the load acts, l is the radius of this zone.
In computations yc = 20 µm, l = 10 µm. At the upper and lower boundaries the
periodicity conditions are given. The size of rectangular computational domain
is 100 µm×40 µm, the dimension of a finite-difference grid is 2560×1024 meshes.
Computations were performed for the liquid crystal 5CB with the parameters:
ρ = 1022 kg/m3 , j = 1.33 · 10−7 kg/m, α = 0.161 GPa, γ = 1 mN, η = 100 Pa · s.
Fig. 5. The action of three Π-shaped impulses of tangential stress at the lateral boun-
daries: level curves of tangential stress at different instants of time.
Figure 6 shows the results of computations for the problem on the action of
a concentrated tangential stress q = q̄ δ(x) δ(t) at the upper boundary. So, at
this boundary: q = q̄, if |x − xc | 6 ∆x (xc = 7 µm is the point of loading)
and t 6 ∆t, and q = 0, otherwise; ω,y = 0. The lower boundary is fixed:
q = 0, ω = 0. At the left and right boundaries the periodicity conditions are
given. The size of computational domain is ten times less than in the previous
problem: 10 µm×4 µm. The dimension of a grid is the same: 2560×1024 meshes.
Computations were performed for the crystal 5CB with parameters as in the
previous case, but j = 1.33 · 10−10 kg/m, γ = 10 µN. In Fig. 6 a the level curves
of tangential stress, and in Fig. 6 b the level curves of angular velocity at different
instants of time are represented.
484
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
a
b
Fig. 6. Lamb’s problem for the action of concentrated tangential stress at the upper
boundary: level curves of tangential stress (a) and angular velocity (b) at different
instants of time.
Figure 7 demonstrates numerical results for the problem on periodic action
of a tangential stress on the part of upper boundary. The boundary conditions
at the upper boundary: q = q̄ sin(2 π ν t), if |x − xc | 6 l, and q = 0, if |x − xc | > l
(xc = 5 µm, l = 2.5 µm); ω,y = 0. As before, the lower boundary is fixed, at the
left and right boundaries the periodicity conditions are given. The frequency ν
is equal to ν∗ = 350 MHz. Computations were performed with parameters as in
the previous case. The fields of angular velocity are shown in Fig. 7.
Fig. 7. Periodic action of tangential stress at the upper boundary: fields of angular
velocity at different instants of time.
485
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
8 Conclusions
Within the framework of acoustic approximation of the model of a liquid crystal,
the system of second-order equations for tangential stress and angular veloci-
ty taking into account the couple-stress interactions was derived. The parallel
computational algorithm for numerical solution of this system was developed.
A comparison of exact and numerical solutions for one-dimensional problem was
fulfilled, and good correspondence of results was obtained. Some computations
were performed by means of the parallel program using the CUDA technology
for computer systems with graphic accelerators.
Acknowledgments. This work was partially supported by the Complex Fun-
damental Research Program no. II.2P “Integration and Development”of SB RAS
(project no. 0356–2016–0728), the Russian Foundation for Basic Research (grants
no. 14–01–00130, 16–31–00078).
References
1. Blinov, L.M.: Structure and Properties of Liquid Crystals. Springer, Heidelberg –
New York – Dordrecht – London (2011)
2. Ericksen, J.L.: Conservation Laws for Liquid Crystals. Trans. Soc. Rheol. 5, 23–34
(1961)
3. Leslie, F.M.: Some Constitutive Equations for Liquid Crystals. Arch. Ration. Mech.
Anal. 28, 265–283 (1968)
4. Aero, E.L., Bulygin, A.N.: The Equations of Motion of Nematic Liquid Crystals.
J. Appl. Math. Mech. 35 (5), 879–891 (1971) [in Russian]
5. Sadovskii, V.M., Sadovskaya, O.V.: On the Acoustic Approximation of Thermome-
chanical Description of a Liquid Crystal. Phys. Mesomech. 16 (4), 312–318 (2013)
6. Sadovskii, V.M.: Equations of the Dynamics of a Liquid Crystal under the Influence
of Weak Mechanical and Thermal Perturbations. AIP Conf. Proc. 1629, 311–318
(2014)
7. Sadovskaya, O.V.: Numerical Simulation of the Dynamics of a Liquid Crystal in
the Case of Plane Strain Using GPUs. AIP Conf. Proc. 1629, 303–310 (2014)
8. Smolekho, I.V.: Parallel Implementation of the Algorithm for Description of Ther-
moelastic Waves in Liquid Crystals. Young Scientist 11 (91), 107–112 (2015)
[in Russian]
9. Samarskiy, A.A., Nikolaev, E.S.: The Methods of Solving Grid Equations. Nauka,
Moscow (1978) [in Russian]
10. Farber, R.: CUDA Application Design and Development. Morgan Kaufmann /
Elsevier, Amsterdam – Boston – Heidelberg – London – New York – Oxford –
Paris – San Diego – San Francisco – Singapore – Sydney – Tokyo (2011)
11. Belyaev, B.A., Drokin, N.A., Shabanov, V.F., Shepov, V.N.: Dielectric Anisotropy
of 5CB Liquid Crystal in a Decimeter Wavelength Range. Phys. Solid State 42 (3),
577–579 (2000)
486