Based on the mathematical model describing the thermomechanical behavior of a liquid crystal, which takes into account the couple-stress interactions, the system of two diferential 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 eficiency of proposed parallel program.
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 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 diferent sections is oriented diferently, the regions with diferent 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 resolution 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 describe 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 dynamics of viscous or inviscid liquid and can rotate relative to a liquid, encountering 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 diferential 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
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 coeficients. Parallel computational algorithm for the solution of this system is represented in [7, 8].
In two-dimensional case the complete system of equations describing the behavior of a liquid crystal under weak acoustic perturbations taking into account the couple stresses is as follows: ρ u,t = −p,x − q,y, j ω,t = 2 q + μx,x + μy,y, φ ,t = ω, ρ v,t = q,x − p,y, p,t = −k( u,x + v,y) + β T,t, q,t = α( v,x − u,y) − 2 α( ω + q/η) , μx,t = γ ω,x,
μy,t = γ ω,y, c T,t = ( ae11 T,x + ae12 T,y) ,x + ( ae12 T,x + ae22 T,y) ,y −
− β T ( u,x + v,y) + 2 q2/η.
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 rotation, η is the viscosity coeficient, c and β are the coeficients of heat capacity and thermal expansion, ae11, ae12 and ae22 are the components of the thermal conductivity tensor: ae11 = ae1 cos2 φ + ae2 sin2 φ , ae12 = (ae1 − ae2) sin φ cos φ , ae22 = ae1 sin2 φ + ae2 cos2 φ , (ae1 and ae2 are the thermal conductivity coefifcients 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 (
Let’s consider how to obtain the system of equations of the second order for
tangential stress and angular velocity. Diferentiating the first equation of (
ρ( 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 diferentiation of corresponding equations of the system (
2 α q,tt + η
2 γ ω,tt − j q,t = j △ω.
α
q,t + 2 α ω,t = ρ △q,
(
x
The fourth-order equation for q can be derived from the subsystem (
) q,t ,
2 γ ω,ttt = j q,tt + j △ω,t, 2 α 4 α ( α γ ) 2 α γ q,tttt + η q,ttt + j q,tt − ρ + j △q,tt − j η Initial data for the corresponding Cauchy problem are as follows: α γ △q,t = − ρ j △2q. q t=0 = q0, q,t t=0 = α( v,0x − u,0y) − 2 α( ω0 + q0 ) = −2 α( ω0 + q0 ) , η η q,tt t=0 = αρ △q0 − 2jα ( 2 q0 + μ0x,x + μ0y,y) − 2ηα q,0t = 4 α[ αη ω0 + ( ηα2 − 1j ) q0], q,ttt t=0 = αρ △q,0t − 2jα
( 2 q,0t + γ △ω0) − 2ηα q,0tt = = 8 α2[( 1 α ) ω0 + 1 ( 1 j − η2 η j + j η − ηα2 ) q0].
1 3
Computational algorithm is developed for numerical solution of the system of two
second-order equations (
Equations of the system (
△t △t = α ( qjn1+1,j2 − 2 qjn1,j2 + qjn1−1,j2 + qjn1,j2+1 − 2 qjn1,j2 + qjn1,j2−1 ) ρ (△x)2 (△y)2 = γ ( ω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, ωjn1+,j12 can be expressed from the second equation: ωjn1+,j12 = 2 ωjn1,j2 − ωjn1−,j12 + △j t ( qjn1+,j12 − qjn1−,j12 )
+
(△x)2
+
γ (△t)2 ( ωjn1+1,j2 − 2 ωjn1,j2 + ωjn1−1,j2 + ωjn1,j2+1 − 2 ωjn1,j2 + ωjn1,j2−1 ) .
j
(△y)2
Substituting (
+ j η △t
+ (△1t)2 ) qjn1+,j12 = (△2t)2 qjn1,j2 +
+ ( αj + η α△t − (△1t)2 ) qjn1−,j12 + 2△αt (
ωjn1−,j12 − ωjn1,j2
)
+
+ α ( qjn1+1,j2 − 2 qjn1,j2 + qjn1−1,j2 + qjn1,j2+1 − 2 qjn1,j2 + qjn1,j2−1 )
ρ (△x)2 (△y)2
−
−
α γ △t ( ωjn1+1,j2 − 2 ωjn1,j2 + ωjn1−1,j2 + ωjn1,j2+1 − 2 ωjn1,j2 + ωjn1,j2−1 ) .
j (△x)2 (△y)2
(
Finite-diference 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.
Under analysis of the stability of the finite-diference 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. According to the Fourier method, let qj1,j2 = λn qˆ ei(j1α1+j2α2), n ωjn1,j2 = λn ωˆ ei(j1α1+j2α2).
Substituting these values into the first equation of the system (
To obtain the characteristic equation, we form the matrix of coeficients under qˆ and ωˆ: (λ(△−t)12)2 + 4ρα λ A
1 λ2 − 1 − j △t α (λ − 1)(λ + 1)
△t (λ2 − 1)2 (△t)2 + 4jγ λ A = 0, where A = sin2(α1/2) (△x)2 + sin2(α2/2) (△y)2
. Introducing the notations a = αρ A (△t)2, b = γj A (△t)2, c = αj (△t)2, 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 diferent values a, b and c: 1) If b = 0, then (1 + c)(λ + 1)2 + 4 λ (a − 1) = 0, λ2 + 2 λ( 1 − 2 1 − a ) + 1 = 0, 1 + c λ1 = λ¯2, | λ1| = | λ2| = 1, ( 1 − 2 11 −+ ac ) 2 −1 6 1 − 2 11 −+ ac 6 1, 2 > 2 11 −+ ac > 0, a 6 1.
Therefore, in this case we obtain the stability condition ∀ α1, α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: αρ (△t)2 6 (
. γj (△t)2 6 (
1 + c | z1| = | z2| = 1, ( The latter condition is satisfied automatically. The condition a + b = 1 means that ( α ρ + γ ) j (△t)2 6 ( sin2(α1/2)
sin2(α2/2) ) −1 + Under such choice of time step, | λ| = 1 for given values α1 and α2.
Thus, in the general case, the following stability condition takes place:
( α
ρ
+
γ )
j
(△t)2 6 (
The described algorithm for numerical solution of the system (
Parallel program has the following structure: 1. Setting the dimensions of finite-diference grid and all the constants used (on the CPU). 2. Description of one-dimensional arrays for tangential stress and angular velocity (on the CPU). 3. Setting the initial data for these variables at the nodes of the finite-diference grid (on the CPU). 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 quantities). 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 procedureskernels 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-diference scheme “cross”; after performing of cores, the barrier synchronization is necessary, to ensure 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 angular 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-diference 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); 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 <<<blocks,threads>>>(it,q,q1,q2,w,w1,w2); cudaThreadSynchronize();
Previously, in solving the problem without taking into account the couplestress interactions, the eficiency of the parallel program for complete system of equations of a model was analyzed [7]. To evaluate the eficiency of parallelization, a large number of computations was performed at diferent grid dimensions. The computation time for parallel program and sequential program was compared. Fig. 1 shows a graph of the dependence of acceleration of the parallel } { } Fig. 1. Acceleration of the program on GPU as compared with CPU program on the dimension N × N of a finite-diference 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
For one-dimensional problem on the action of tangential stress q = qˆ ei(ft−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(ft−ky), ω = ωˆ ei(ft−ky) in the equations of system (
Calculating the determinant of this system, one can find the expression for k±: v uu ρ j f ( k± = t 2 α γ
√ d ± d2 − 4 αρ jγ ( f 2 − 2 i α f η .
Here f is the frequency, k± = k1± + i k2± are the wave numbers.
The characteristic dispersion curves are represented in Figs. 2 and 3: def f pendence of the phase velocity c± = on the frequency ν = and Re k± 2 π 1 dependence of the damping decrement λ± = − Im k± on the frequency ν (for 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 parameters [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 determined from the next equations:
Re qˆ = ek2y( qˆ1 cos(f t − k1y) − qˆ2 sin(f t − k1y)) ,
Re ωˆ = ek2y( ωˆ1 cos(f t − k1y) − ωˆ2 sin(f t − k1y)) .
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-diference grid is 1000 meshes. The relative error is 3 · 10−3 in calculations for k+ and 5 · 10−4 in calculations for k−. 450 400 350 ,μm300 +λ250 200 150
0 2 a1 P ,G0 ω−1 −2 a a a 0 0.5 1 1.5 2 2.5 3 3.5 4 y, μm 7
A series of numerical calculations was carried out on the high-performance computational 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 eficiency 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-diference 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.
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 diferent instants of time are represented. Mathematical and Information Technologies, MIT-2016 | Mathematical modeling a b
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.
Within the framework of acoustic approximation of the model of a liquid crystal, the system of second-order equations for tangential stress and angular velocity 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 Fundamental 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).