A Linear, Decoupled and Energy stable scheme for smectic-A Liquid Crystal Flows

1.3 MB·

by  Xiaofeng Yang

A LINEAR, DECOUPLED AND ENERGY STABLE SCHEME FOR
SMECTIC–A LIQUID CRYSTAL FLOWS
7 XIAOFENGYANG†? ANDALEXBRYLEV‡
1
0
2

n
Abstract. Inthispaper,weconsidernumericalapproximationsforthemodelofsmectic-
a
Aliquidcrystalflows. Themodelequation,thatisderivedfromthevariationalapproach
J
of the de Gennes free energy, is a highly nonlinear system that couples the incompress-
5
2 ible Navier-Stokes equations, and two nonlinear coupled second-order elliptic equations.
Based on some subtle explicit–implicit treatments for nonlinear terms, we develop a un-

] conditionally energy stable, linear and decoupled time marching numerical scheme. We
A
also rigorously prove that the proposed scheme obeys the energy dissipation law at the
N discrete level. Various numerical simulations are presented to demonstrate the accuracy
. andthestabilitythereafter.
h
t
a
m
1. Introduction
[

1 Liquid crystal (LC) is often viewed as the fourth state of the matter besides the gas,
v
liquidandsolid. Itmayflowlikealiquid, butitsmoleculesmaybeorientedinacrystal-like
7
way. There are many different types of liquid-crystal phases, which can be distinguished
8
4 by their different optical properties. Thermotropic LCs can be distinguished into two main
7
different phases: Nematic and Smectic. In Nematic phases, the rod-shaped molecules have
0
. no positional order, but molecules self-align to have a long-range directional order with
1
0 their long axes roughly parallel. Thus, the molecules are free to flow and their center of
7 mass positions are randomly distributed as in a liquid, although they still maintain their
1
long-range directional order. In smectic phases, which are found at lower temperatures, the
:
v
well-definedlayersform, thatcanslideoveroneanotherinamannersimilartothatofsoap.
i
X The smectics are thus positionally ordered along one direction inside the layer. There are
r many different smectic phases, all characterized by different types and degrees of positional
a
andorientationalorder(cf.[3,7,11,12,17,37,39,42–44,44,47–49]andthereferencestherein).
In particular, in Smectic-A phase, molecules are oriented along the normal vector of the
layers, while in Smectic-C phases they are tilted away from the normal vector of the layer.
Key words and phrases. smectic-A, Liquid Crystal, Unconditional Stability, Ginzburg-Landau, Decoupled,
Linear.
†? Corresponding author, Department of Mathematics, University of South Carolina, Email:
[email protected].
‡ DepartmentofMathematics,UniversityofSouthCarolina,Email: [email protected].
1
2 XIAOFENGYANGANDALEXBRYLEV
Themathematicalmodelofliquidcrystalscanoftenbederivedfromanenergy-basedvari-
ational formalism (energetic variational approach), leading to well-posed nonlinearly cou-
pled systems that satisfy thermodynamics-consistent energy dissipation laws. This makes
it possible to carry out mathematical analysis and design numerical schemes which satisfy
a corresponding energy dissipation law at the discrete level. For smectic-A phase liquid
crystals, indeGennes’pioneeringwork[7], thephenomenonlogicalfreeenergyofsmectic–A
phase is presented by coupling two order parameters which characterize the average direc-
tion of molecular alignment, as well as the layer structure, respectively. In [4], the authors
modified the de Gennes’ model by adding a second order gradient term for the smectic
order parameter to investigate the nematic to smectic–A or smectic–C phase transition,
and to predict the twist grain boundary phase in chiral smectic liquid crystals. In [18],
the authors used the de Gennes energy to study smectic–A liquid crystals to simulate the
chevron (zigzag) pattern formed in the presence of an applied magnetic field. In [8], the
authors derived the hydrodynamics coupled model for smectic–A phase by assuming that
thedirectorfieldisstrictlyequaltothegradientofthelayer, thusthefreeenergyisreduced
to one order parameter.
From the numerical point of view, it is specifically desired to design numerical schemes
thatcouldpreservethethermo-dynamicallyconsistentdissipationlaw(energy-stable)atthe
discretelevel,sincethepreservationofsuchlawsiscriticalfornumericalmethodstocapture
the correct long time dynamics. The noncompliance of energy dissipation laws may lead to
spurious numerical solutions if the grid and time step sizes are not carefully controlled. To
the best of the author’s knowledge, although a variety of the smectic liquid crystal models
had been developed, we notice that the successful attempts in designing efficient energy
stable schemes are very scarce due to the complex nonlinearities. For instances, in [14], the
authors present a temporal second order numerical scheme to solve the model of [8]. The
schemeisenergystable,however,itisnonlinearthustheimplementationiscomplicatedand
the computational cost might be high. In [18], the authors developed a temporal first order
scheme. However, it does not follow the energy dissipation law even though the schemes
are linear and decoupled.
Therefore,themainpurposeofthispaperistoconstructtheefficientschemestosolvethe
deGennestypesmectic–Aliquidcrystalmodel(cf.[7,18]). Wefirstcouplethehydrodynam-
icstotheoriginaldeGennesfreeenergyandderivethewholemodelbasedonthevariational
approach and the Ficks’ law. To solve the model, the main difficulties roughly include (i)
the coupling between the velocity and director field/layer function through the convection
terms and nonlinear stresses; (ii) the coupling of the velocity and pressure through the in-
compressibility constraint; (iii) the nonlinear coupling between the director field and layer
function. We develop a time discretization scheme which (a) is unconditionally stable; (b)
satisfies a discrete energy law; and (c) leads to linear, decoupled equations to solve at each
time step. This is by no means an easy task due to many highly nonlinear terms and the
couplings existed in the model.
NUMERICAL APPROXIMATIONS FOR SMECTIC–A LIQUID CRYSTAL FLOWS 3
The rest of the paper is organized as follows. In Section 2, we present the whole model
and present the PDE energy law. In Section 3, we develop the numerical scheme and prove
theunconditionalstability. InSection4,wepresentsomenumericalexperimentstovalidate
the proposed scheme. Finally, some concluding remarks are presented in Section 5.
2. The smectic-A liquid crystal fluid flow model and its energy law
The de Gennes free energy of smectic A liquid crystal is described by a unit vector
(director field) d and a complex order parameter ψ, to represent the average direction of
molecular alignment and the layer structure, respectively. The smectic order parameter is
written as
(2.1) ψ(x)=ρ(x)eiqω(x),
whereω(x)istheorderparametertodescribethelayerstructuresothat ωisperpendicular
∇
tothelayer. Thesmecticlayerdensityρ(x)isthemassdensityofthelayers. ThedeGennes
free energy reads as follows,
g r
(2.2)E(ψ,d)= C ψ iqdψ 2+K d2+ (ψ 2 )2 χ H2(d h)2 dx,
a
|∇ − | |∇ | 2 | | − g − ·
ZΩ(cid:16) (cid:17)
where the order parameters C,K,g,r are all fixed positive constants, h is the unit vector
representing the direction of magnetic field and H2 is the strength of the applied field.
Ω=( L,L)2 ( d,d).
− × −
Nowweconsiderthesimplecasebyassumingthedensityρ(x)=r/g [18],thentheenergy
becomes
(2.3) E(ψ,d)= Cq2 ω d2+K d2 χ H2(d h)2 dx.
a
|∇ − | |∇ | − ·
ZΩ(cid:16) (cid:17)
ω(x)
Let φ(x)= , thus the normalized energy becomes
d
1 φ d2 d2 τ
(2.4) E(φ,d)=λ |∇ − | +η|∇ | (d h)2 dx,
ZΩ˜ (cid:16)η 2 2 − 2 · (cid:17)
where
x L
x˜ = , Ω˜ =(0,2‘)2 (0,2), ‘= ,
d × d
(2.5)
γ K 2dK χ H2d2η
a
η = ,γ = , λ= , τ = .
d sCq2 η K
The dimensionless parameter η is in fact the ratio of the layer thickness to the sample
thickness and thus η 1.
(cid:28)
Toreleasetheunitvectorconstraintof d =1,anonlinearpotentialG(d)= 1 (d2 1)2
| | 4(cid:15)2 | | −
whichisaGinzburg-Landautypepenaltyterm,isaddedintothefreeenergytoapproximate
the unit length constraint of d [20,21] , where (cid:15) 1 is a penalization parameter. Thus the
(cid:28)
4 XIAOFENGYANGANDALEXBRYLEV
modified total free energy with the hydrodynamics reads as follows
1 1 φ d2 d2 τ
(2.6)E (u,φ,d)= u2dx+λ |∇ − | +η|∇ | +G(d) (d h)2 dx,
tot
ZΩ˜ 2| | ZΩ˜ (cid:16)η 2 2 − 2 · (cid:17)
where u is the fluid velocity field.
Assuming a generalized Fick’s law that the mass flux be proportional to the gradient of
the chemical potential [2,10,22,24], we can derive the following system:
δE
(2.7) φ +u φ= M µ , µ = ,
t 1 φ φ
·∇ − δφ
δE
(2.8) d +u d= M µ , µ = ,
t ·∇ − 2 d d δd
(2.9) u +u u σ + p µ φ µ d=0,
t d φ d
·∇ −∇· ∇ − ∇ − ∇
(2.10) u=0,
∇·
where p is the pressure, σ is the Caughy stress tensor, ν is the viscosity, M ,M are the
d 1 2
relaxation order parameters. The variational derivatives µ and µ are
φ d
1
(2.11) µ =λ ( ∆φ+ d),
φ
η − ∇·
1
(2.12) µ =λ η∆d+g(d)+ ( φ+d) τ(d h)h ,
d − η −∇ − ·
(cid:16) (cid:17)
where g(d)= 1d(d2 1).
(cid:15)2 | | −
Following the work in [8], the Cauchy stress tensor σ reads as follows,
d
(2.13) σ =µ (dTD(u)d)d d+µ D(u)+µ (D(u)d d+d D(u)d),
d 1 4 5
⊗ ⊗ ⊗
where µ >0 and D(u)= 1(( u)+( u)T) is the strain tensor. In this paper, we assume
i 2 ∇ ∇
µ ,µ are negligible comparing to µ , thus the stress tensor is simplifed to
1 5 4
(2.14) σ =µ D(u).
d 4
For simplicty, we assume the boundary conditions as follows.
(2.15) u =0,∂ φ =0,∂ d =0,
∂Ω n ∂Ω n ∂Ω
| | |
where n is the outward normal of the boudary.
To obtain the dissipation law of the system (2.7)-(2.9), we take the L2 inner product of
(2.7) with δE, (2.8) with δE, and (2.9) with u, perform the integration by parts, and add
δφ δd
all equalities together, we obtain
d δE δE
(2.16) E (u,φ,d)= µ D(u)2+M 2+M 2 dx 0.
tot 4 1 2
dt − | | |δφ| |δd| ≤
ZΩ(cid:16) (cid:17)
NUMERICAL APPROXIMATIONS FOR SMECTIC–A LIQUID CRYSTAL FLOWS 5
3. Numerical scheme
The emphasis of our algorithm development is placed on designing numerical schemes
that are not only easy-to-implement, but also satisfy a discrete energy dissipation law. We
will design schemes that in particular can overcome the following difficulties, namely,
the coupling of the velocity and pressure through the incompressibility condition;
•
the stiffness in the director equation associated with the penalty parameter (cid:15);
•
thenonlinearcouplingsamongthefluidequation,thelayerequationandthedirector
•
equation.
We construct an energy stable scheme based on a stabilization approach [30,31]. To this
end, we shall assume that G(d) satisfies the following conditions, i.e.,
(3.1) HG(x) L, x R3.
| |≤ ∀ ∈
where (H (x)) = ∂2G ,i,j = 1,2,3 is the Hessian matrix of G(x). One immediately
G i,j ∂xi∂xj
notes that this condition is not satisfied by this double-well potential G(x). However,
it is a common practice that one truncates this fourth order polynomial G to quadratic
growth outside of an interval [ M,M] without affecting the solution if the maximum norm
−
of the initial condition d is bounded by M. Therefore, one can (cf. [6,19,30]) consider
0
| |
the truncated double-well potential G˜(d) by modifying this function outside a ball in x :
{
d(x) 1 R3 of radius 1 as follows.
| |≤ }∈
1
(d2 1)2, d 1,
(3.2) G˜(d)= 4(cid:15)2 | | − | |≤
 1
 (d 1)2, d >1.
4(cid:15)2 | |− | |
Hence, there exists a postive constant LG such that
(3.3) max H (x) L.
x R3| G˜ |≤
∈
For convenience, we consider the problem formulated with the substitute G˜, but omit the
˜in the notation.
Whenderivingtheenergylaw (2.16), wenoticethatthenonlineartermsin δEtot and δEtot
δφ δd
involve second order derivatives, and it is not convenient to use them as test functions in
numerical approximations, making it difficult to prove the discrete energy dissipation law.
To overcome this difficulty, we first reformulate the system (2.7)-(2.10) in an alternative
form which is convenient for numerical approximations. The system reads as follows,
(3.4) φ +u φ= M µ ,
t 1 φ
·∇ −
(3.5) d +u d= M µ ,
t ·∇ − 2 d
φ˙ d˙
(3.6) u +u u σ + p+ φ+ d=0,
t d
·∇ −∇· ∇ M ∇ M ∇
1 2
(3.7) u=0,
∇·
6 XIAOFENGYANGANDALEXBRYLEV
where φ˙ = φ +u φ and d˙ = d +u d. To obtain the dissipation law of the system
t t
·∇ ·∇
(3.4)-(3.7), we take the L2 inner product of (3.4) with φ , (3.5) with d , and (3.6) with u,
t t
perform the integration by parts, and add all equalities together. We obtain
d 1 1
(3.8) E = µ u2+ φ˙ 2+ d˙ 2 dx 0.
tot 4
dt − |∇ | M | | M | | ≤
ZΩ(cid:16) 1 2 (cid:17)
We now fix some notations. For scalar function u,v and vector function u=(u ,u ,u )
1 2 3
and v =(v ,v ,v ), we denote the L2 inner product as follows.
1 2 3
(3.9) (u,v)= uvdx, u 2 =(u,u), (u,v)= uvTdx, u 2 =(u,u).
k k k k
ZΩ ZΩ
Now, we are ready to present our energy stable scheme that reads as follows.
Given the initial conditions φ0, d0, u0 and p0 = 0, having computed φn, dn, un and pn
for n>0, we compute (φn+1,dn+1,u˜n+1,un+1,pn+1) by
Step 1:
1 λ
φ˙n+1 = (∆φn+1 dn),
M η −∇·
(3.10)  1
∂φn+1
 =0,
∂Ω
∂n |
with 

φn+1 φn φ˙n+1
(3.11) φ˙n+1 = − +(un )φn, un =un δt φn.
δt ? ·∇ ? − M ∇
1
Step 2:
S(dn+1 dn)+ 1 d˙n+1
− M
2

1
(3.12) =λ η∆dn+1 g(dn)+ ( φn+1 dn+1)+τ(dn h) h ,
 − η ∇ − · ·

 (cid:16) (cid:17)
∂dn+1
=0,
 ∂n |∂Ω


with

(3.13) d˙n+1 = dn+1−dn +(un )dn, un =un δtd˙n+1 dn.
δt ??·∇ ?? ? − M ∇
2
Step 3:
u˜n+1 un φ˙n+1 d˙n+1
− +(un )u˜n+1 µ ∆u˜n+1+ pn+ φn+ dn =0,
4
(3.14) δt ·∇ − ∇ M ∇ M ∇
 1 2
u˜n+1 =0,
∂Ω
|


NUMERICAL APPROXIMATIONS FOR SMECTIC–A LIQUID CRYSTAL FLOWS 7
Step 4:
un+1 u˜n+1
− + (pn+1 pn)=0,
δt ∇ −
(3.15)  un+1 =0,
∇·

n un+1 =0.
∂Ω
· |

In the above, S is a stabilizing parameter to be determined.

We have the following remarks in order:
Apressure-correctionscheme[13]isusedtodecouplethecomputationofthepressure
•
from that of the velocity.
The nonlinear term g(d) mainly takes the form like 1d(d2 1), so the explicit
• (cid:15)2 | | −
treatment of this term usually leads to a severe restriction on the time step δt when
(cid:15) 1. Thus we introduce in (3.10) a linear “stabilizing” term to improve the
(cid:28)
stability while preserving the simplicity. It allows us to treat the nonlinear term
explicitly without suffering from any time step constraint [5,23,25,28–30,32–36,38,
40,41,45–47]. Note that this stabilizing term introduces an extra consistent error of
order O(δt) in a small region near the interface, but this error is of the same order
as the error introduced by treating it explicitly, so the overall truncation error is
essentially of the same order with or without the stabilizing term. It is remarkable
that the truncation error of the stabilizing approach is exactly the same as the
convex splitting method [9].
Inspired by [1,26,31], which deal with a phase-field model of three-phase viscous
•
fluidsorcomplexfluids, weintroducetwonew, explicit, convectivevelocitiesun and
?
un in the phase equations. un and un can be computed directly from (3.11) and
?? ? ??
(3.13), i.e.,
δt 1 1
(3.16) un = I + ( φn)T φn − un (φn+1 φn) φn ,
? M ∇ ∇ − M − ∇
1 1
(cid:16) (cid:17) (cid:16) (cid:17)
δt 1 1
(3.17) un = I + ( dn)T dn − un (dn+1 dn) dn .
?? M ∇ ∇ ? − M − ∇
2 2
(cid:16) (cid:17) (cid:16) (cid:17)
It is easy to get det(I + c( φ)T φ) = 1 + c φ φ, thus the above matrix is
∇ ∇ ∇ · ∇
invertible.
Thescheme(3.10)-(3.15)isatotallydecoupled,linearscheme. Indeed,(3.10),(3.12)
•
and (3.14) are respectively (decoupled) linear elliptic equations for φn+1, dn+1 and
u˜n+1,and(3.15)canbereformulatedasaPoissonequationforpn+1 pn. Therefore,
−
at each time step, one only needs to solve a sequence of decoupled elliptic equations
that can be solved very efficiently.
As we shall show below, the above scheme is unconditionally energy stable.
•
8 XIAOFENGYANGANDALEXBRYLEV
Theorem 3.1. Under the condition (3.3), and S λL, the scheme (3.10)-(3.15) admits a
≥ 2
unique solution satisfying the following discrete energy dissipation law:
δt2 φ˙n+1 2 d˙n+1 2 δt2
En+1+ pn+1 2+ νδt u˜n+1 2+δt(| | + | | ) En+ pn 2,
2 k∇ k k∇ k M M ≤ 2 k∇ k
(cid:26) 1 2 (cid:27)
where
1 dn 2 1 dn φn 2 τ
(3.18)En = un 2+λ ηk∇ k +(G(dn),1)+ k −∇ k dn h 2 .
2k k 2 η 2 − 2k · k
(cid:16) (cid:17)
Proof. Fromthedefinitionofun andun in(3.11)and(3.13),wecanrewritethemomentum
? ??
equation (3.14) as follows
u˜n+1 un
(3.19) − ?? +(un )u˜n+1 µ ∆u˜n+1+ pn =0.
4
δt ·∇ − ∇
By taking the inner product of (3.19) with 2δtu˜n+1, and using the identity
(3.20) (a b,2a)= a2 b2+ a b2,
− | | −| | | − |
we obtain
(3.21)u˜n+1 2 un 2+ u˜n+1 un 2+2µ δt u˜n+1 2+2δt( pn,u˜n+1)=0.
k k −k ??k k − ??k 4 k∇ k ∇
To deal with the pressure term, we take the inner product of (3.15) with 2δt2 pn to derive
∇
(3.22) δt2( pn+1 2 pn 2 pn+1 pn 2)=2δt(u˜n+1, pn).
k∇ k −k∇ k −k∇ −∇ k ∇
By taking the inner product of (3.15) with un+1, we obtain
(3.23) un+1 2+ un+1 u˜n+1 2 = u˜n+1 2.
k k k − k k k
We also derive from (3.15) directly that
(3.24) δt2 pn+1 pn 2 = u˜n+1 un+1 2.
k∇ −∇ k k − k
Combining all identities above, we obtain
un+1 2 un 2+ u˜n+1 un 2
(3.25) k k −k ??k k − ??k
+δt2( pn+1 2 pn 2)+2νδt u˜n+1 2 =0.
k∇ k −k∇ k k∇ k
Next, we derive from (3.11) and (3.13) that
un un φ˙n+1
(3.26) ? − = φn,
δt − M ∇
1
un un d˙n+1
(3.27) ??− ? = dn.
δt − M ∇
2
NUMERICAL APPROXIMATIONS FOR SMECTIC–A LIQUID CRYSTAL FLOWS 9
By taking the inner product of (3.26) with 2δtun, of (3.27) with 2δtun , we obtain
? ??
φ˙n+1
(3.28) un 2 un 2+ un u 2 = 2δt( φn,un),
k ?k −k k k ? − nk − M ∇ ?
1
d˙n+1
(3.29) un 2 un 2+ un un 2 = 2δt( dn,un ).
k ??k −k ?k k ??− ?k − M ∇ ??
2
Then, by taking the inner product of (3.10) with 2(φn+1 φn), we obtain
−
φ˙n+1 2 φ˙n+1 λ
2δtk k 2δt( ,(un )φn)+ φn+1 2 φn 2+ φn+1 φn 2
M − M ? ·∇ η k∇ k −k∇ k k∇ −∇ k
(3.30) 1 1
(cid:16) (cid:17)
2λ
+ dn,φn+1 φn =0.
η ∇· −
(cid:16) (cid:17)
By taking the inner product of (3.12) with 2(dn+1 dn), we arrive at
−
d˙n+1 2 d˙n+1
2S dn+1 dn 2+2δtk k 2δt( ,(un )dn)
k − k M − M ??·∇
2 2
dn+1 2 dn 2 dn+1 dn 2
+2λη k∇ k k∇ k + k∇ −∇ k
2 − 2 2
(3.31) + λ (cid:16)dn+1 2 dn 2+ dn+1 dn 2 (cid:17)
η k k −k k k − k
(cid:0) 2λ (cid:1)
+2λ(g(dn),dn+1 dn) φn+1,dn+1 dn
− − η ∇ −
(cid:16) (cid:17)
2λτ (dn h)h,dn+1 dn =0.
− · −
(cid:16) (cid:17)
10 XIAOFENGYANGANDALEXBRYLEV
Combining (3.25), (3.28), (3.29), (3.30), and (3.31), we arrive at
un+1 2 un 2+ u˜n+1 un 2+ un un 2+ un un 2
k k −k k k − ??k k ??− ?k k ? − k
+δt2( pn+1 2 pn 2)
k∇ k −k∇ k
φ˙n+1 2 d˙n+1 2
+2νδt u˜n+1 2+2δtk k +2δtk k
k∇ k M M
1 2
dn+1 2 dn 2 dn+1 dn 2
+2λη k∇ k k∇ k + k∇ −∇ k
2 − 2 2
λ (cid:16) (cid:17)
+ ( dn+1 2 dn 2+ dn+1 dn 2)
η k k −k k k − k
(3.32) λ
+ ( φn+1 2 φn 2+ φn+1 φn 2)
η k∇ k −k∇ k k∇ −∇ k
+2S dn+1 dn 2
k − k
+2λ(g(dn),dn+1 dn)(:Term A)
−
2λ 2λ
+ ( dn,φn+1 φn) ( φn+1,dn+1 dn)(:Term B)
η ∇· − − η ∇ −
2λτ (dn h)h,dn+1 dn (:Term C)
− · −
=0. (cid:16) (cid:17)
We deal with the terms A,B,C as follows.
For Term A, we apply the Taylor expansions to obtain
g (ξ)
(3.33) A=2λ(G(dn+1) G(dn),1) 2λ 10 , dn+1 dn 2 .
− − 2 | − |
(cid:16) (cid:17)
For Term B, we have
λ
B = 2 (dn, φn+1 φn)+(dn+1 dn, φn+1)
− η ∇ −∇ − ∇
(3.34) (cid:16) (cid:17)
λ
= 2 (dn+1, φn+1) (dn, φn) .
− η ∇ − ∇
(cid:16) (cid:17)
For Term C, we have
C = 2λτ (dn h),(dn+1 h) (dn h)
− · · − ·
(3.35)
(cid:16) (cid:17)
= λτ dn+1 h 2 dn h 2 (dn+1 dn) h 2 .
− k · k −k · k −k − · k
(cid:16) (cid:17)