Finite element approximation of nematic liquid crystal flows using a saddle-point structure

FINITE ELEMENT APPROXIMATION OF NEMATIC LIQUID CRYSTAL FLOWS USING A SADDLE-POINT STRUCTURE SANTIAGO BADIA∗, FRANCISCO GUILL´ EN-GONZ´ ALEZ †,AND JUAN VICENTE GUTI´ ERREZ-SANTACREU ‡ Abstract. In this work, we propose finite element schemes for the numerical approximation of nematic liquid crystal flows, based on a saddle-point formulation of the director vector sub-problem. It introduces a Lagrange multiplier that allows to enforce the sphere condition. In this setting, we can consider the limit problem (without penalty) and the penalized problem (using a Ginzburg-Landau penalty function) in a unified way. Further, the resulting schemes have an stable behavior with respect to the value of the penalty parameter, a key difference with respect to the existing schemes. Two different methods have been considered for the time integration. First, we have considered an implicit algorithm that is unconditionally stable and energy preserving. The linearization of the problem at every time step value can be performed using a quasi-Newton method that allows to decouple fluid velocity and director vector computations for every tangent problem. Then, we have designed a linear semi-implicit algorithm (i.e. it does not involve nonlinear iterations) and proved that it is unconditionally stable, verifying a discrete energy inequality. Finally, some numerical simulations are provided. Key words. Nematic liquid crystals, finite element methods, saddle-point problems, Ericksen-Leslie problem, Ginzburg-Landau problem 1. Introduction. Liquid crystals are materials that exhibit intermediate transitions from liquid to solid made of anisotropic macro-molecules of similar size. Typically, liquid crystal molecules are represented as rods. The nematic phase is considered the simplest liquid crystal phase where the elongatedmolecules tend to be locally parallel to some preferentialdirection. However the molecular centers of gravity are allowed to flow freely as in an isotropic fluid, i.e. without a positional order. This uniaxial orientational order is typically modelled by a unit vector called the director vector d. The first phenomenologicaltheory describing statical configurations of a nematic liquid crystal was proposed by Oseen [27] and Frank [15]. They suggested that the director field corresponds to a minimum of the so-called Oseen-Frankfree- energy functional, which in the most elementary form is the Dirichlet energy (1.1) E(d) = KZΩ |∇d|2dx subject to the sphere condition; Kis an elastic constant. It is known that orientational orders affect all the macroscopic properties of the fluid velocity, introducing an anisotropic stress tensor in the linear momentum equations. The ∗International Center for Numerical Methods in Engineering (CIMNE), Universitat Polit`ecnica de Catalunya, Jordi Girona 1-3, Edifici C1, 08034 Barcelona, Spain. [email protected]. †Dpto. de Ecuaciones Diferenciales y An´alisis Num´erico, Universidad Sevilla. Avda. Reina Mercedes, s/n. 41012 Sevilla, Spain. [email protected] ‡Dpto. de Matem´atica Aplicada I, Universidad Sevilla. Avda. Reina Mercedes, s/n. 41012 Sevilla, Spain. [email protected] 1 2 continuum theory of nematic liquid crystals was formulated by Ericksen [11, 12] and Leslie [21, 20], containing the Oseen-Frank elastic energy. Our interest is to construct numerical approximations for the motion of a nematic liquid crystal governed by the simplification of the Ericksen-Leslie equations proposed by Lin in [22]. This problem has numerically been treated using the Ginzburg-Landau penalty problem in order to enforce the sphere constraint. In Section 2 we formulate the problem in the saddle-point framework. Such a formulation allows us to enforce the sphere condition with or without penalty in a single setting. One benefit of this approach is that an energy estimate is obtainedfor both cases. In Section 3 we present a semi-discrete scheme based on low-orderfi- nite elements for approximating all the unknowns. This scheme is unconditionally stable and its solution satisfies a discrete energy estimate. In Section 4, three time-stepping schemes are considered. The first two schemes are nonlinear, with a backward-differencingand mid-point discretization, respectively. Different linearizations for these schemes are studied in Section 5. The third scheme is linear, implicit with respect to the linear term and semi-implicit with respect to the nonlinear term. The three schemes are again unconditionally stable. In Section 6, we test our numerical algorithms with an initial smooth condition and a initial condition with two defect points. Finally, we compare the numerical approximations for a test with analytical solution. 2. Problem statement. A micro-macroscopic continuum theory has been developed for the modeling of nematic liquid crystal flows (see [10]), that characterizes this physical phenomenon in terms of the (microscopic) molecular orientation and the (macroscopic) velocity-pressure variables. The simplified Ericksen-Leslie system consists of a set of partial differential equations that reads as follows: find d,u, and ˜psuch that (2.1) ∂td+ (u· ∇)d−γ∆d−γ|∇d|2d=0, |d|= 1, ∂tu+ (u· ∇)u−ν∆u+∇ep+λ∇ · ((∇d)t∇d) = g, ∇ · u= 0, in (x, t)∈Ω×(0, T ), where Ω⊂R3is the spatial bounded domain filled by the liquid crystal, and [0, T ]the time interval. The physical constants are the fluid viscosity ν > 0, an elasticity constant λ > 0and a relaxation time γ > 0. The unknown d(x, t)∈R3is the director vector that determines the orientation of the molecules, u(x, t)∈R3is the velocity of the flow and ep(x, t)∈Rits pressure. The data g(x, t)∈R3is a force term. The gradient operator is defined as ∇x=∂jxiand (∇x)tdenotes its transpose. In the following, we will consider that the boundary conditions ∂nd= 0 and u= 0 are satisfied a.e. on the boundary ∂Ωif we do not specify the contrary; ∂nd=∇d·nis the normal derivative where nis the outward normal vector to the boundary. Initial boundary conditions d(x,0) = d0(with |d0|= 1 a.e. in Ω), and u(x,0) = u0(with ∇ · u0= 0 in Ωand u0= 0 on ∂Ω). 3 In fact, system (2.1) is a simplification of the classical Ericksen-Leslie theory of liquid crystals obtained after assuming that some physical elastic constants are equal (see [10]). In general, this assumption is not true, but the mathematical nature of the system does not change, and the complications related to its numerical approachare still present in the simplified problem. For this reason, system (2.2) has been subject of many mathematical analyses (see [17, 25, 26, 4]). The saddle-point formulation for (2.1) consists in finding d,u,q, and ˜psuch that ∂td+ (u· ∇)d+γ(−∆d+qd) = 0, d·d= 1,(2.2) ∂tu+ (u· ∇)u−ν∆u+λ∇ · ((∇d)t∇d) + ∇˜p=g, ∇ · u= 0, where q(x, t)∈Ris the Lagrange multiplier used to enforce the restriction |d|2= 1 a.e. in Ω×(0, T )(| · | denotes the Euclidean norm). It is interesting to observe that problem (2.2) has only quadratic nonlinear terms whereas problem (2.1) involves cubic nonlinear terms. In top of the open questions related to the Navier-Stokes equations, the nonconvex constraint over dmakes the theoretical analysis of the previous problem very difficult to approach. So, a penalized version is usually considered, in which the constraint |d|2= 1 is weakly enforced by adding the Ginzburg-Landau (GL) function fǫ(d), for 0< ǫ ≪1, where fǫ(d) := 1 ǫ2(d·d−1)d. We will also make use of the potential function Fǫ(d) := 1 4ǫ2(d·d−1)2; note that ∇dFǫ(d) = fǫ(d). The GL penalized problem then reads ∂td+ (u· ∇)d+γ(−∆d+fǫ(d)) = 0, ∂tu+ (u· ∇)u−ν∆u+λ∇ · ((∇d)t∇d) + ∇˜p=g,(2.3) ∇ · u= 0, supplemented with the respectiveinitial and boundaryconditions. In fact, it is straightforward to note that the GL penalized problem (2.3) can also be casted in a saddle-point form as follows: ∂td+ (u· ∇)d+γ(−∆d+qd) = 0, d·d−ǫ2q= 1,(2.4) ∂tu+ (u· ∇)u−ν∆u+λ∇ · ((∇d)t∇d) + ∇˜p=g, ∇ · u= 0, where q=1 ǫ2(d·d−1), hence qd=fǫ(d). The clear advantage of (2.4) with respect to (2.3) is the fact that it formally allows ǫ= 0, i.e. it includes the limit and penalized problem in a unified formulation. This unified approach also permits to make connections between 4 existing methods that seemed essentially different. Further, the saddle-point approach gives a clue about how to deal with the mathematical analysis of the limit problem. The stability of the multiplier qin the limit case can only be attained via an inf-sup conditions ` a la Babuska– Brezzi, with a formidable complication: the inf-sup condition is nonlinear. For the subsequent numerical analysis, we will consider a reformulation of the coupling term in the fluid momentum equation. After some manipulation, the coupling term can be re-written as: ∇·((∇d)t∇d) = ∂j(∂jdk∂idk) = ∂2 jdk∂idk+1 2∂i(|∂jdk|2) = (∇d)t∆d+1 2∇(|∇d|2). Using the second equation in (2.4), we have: q(∇d)td=1 2q∇(|d|2) = 1 2q∇(ǫ2q) = ǫ2 4∇q2. Using the d-system in (2.4), we get: ∇ · ((∇d)t∇d) = (∇d)t(∆d−qd) + 1 2∇(|∇d|2+ǫ2 2q2) =1 γ(∇d)t(∂td+ (u· ∇)d) + 1 2∇(|∇d|2+ǫ2 2q2). Note that for (2.3) one obtains q(∇d)td=q1 2∇(|d|2) = 0. Then, we have the above equality for ǫ= 0. We can absorb the second term in the pressure gradient by using the modified pressure p= ˜p+λ 2|∇d|2+ǫ2 4q2, leading to the following version of the fluid momentum equation: (2.5) ∂tu+ (u· ∇)u−ν∆u+λ γ(∇d)t(∂td+ (u· ∇)d) + ∇p=g. In fact, these manipulations are not new for the GL penalized problem (see e.g. [17, 4]). This reformulation will allow us to obtain finite element approximations with an energy estimate that mimics the one at the continuous level. Since the aim of this work is to consider a Faedo-Galerkin approximationof system (2.3) or (2.4) based on finite element spaces, we consider the problem in a weak sense as follows: find (d(t), q(t),u(t), p(t)) ∈W1,3(Ω) ∩L∞(Ω)×H−1 ǫ(Ω)×H1 0(Ω)×L2 0(Ω) such that (∂td,¯ d) + γ(∇d,∇¯ d) + c(u,d,¯ d) + γbd(q, d,¯ d) = 0,(2.6a) bd(¯q, d,d)−ǫ2(q, ¯q) = h1,¯qi,(2.6b) (∂tu,¯ u) + ν(∇u,∇¯ u) + c(u,u,¯ u) +λ γc(¯ u,d, ∂td+ (u· ∇)d)−bu(p, ¯ u) = hg,¯ ui,(2.6c) bu(¯p, u) = 0,(2.6d) 5 hold for any (¯ d,¯q, ¯ u,¯p)∈H1(Ω)×H−1 ǫ(Ω)×H1 0(Ω)×L2 0(Ω) a.e. in t∈(0, T ). The forms associated to this problem are defined as bd(q, d,¯ d) = (q, d·¯ d), bu(p, ¯ u) = (p, ∇ · ¯ u), c(u,v,¯ u) = h(u· ∇)v,¯ ui. Note that the elastic tensor effect in the u-system is denoted by the same form c(·,·,·). Hereafter, (·,·)denotes the inner product in L2(Ω), with k · k the associated norm, and h·,·i the duality product between (H1(Ω))′and H1(Ω). Since the Lagrange multiplier qwill lose regularity from ǫ > 0to ǫ= 0, we consider the Banach space H−1 ǫ(Ω); for any ǫ > 0it is the Hilbert space L2(Ω) but endowed with the norm ǫk · k +k · k−1, whereas for ǫ= 0 it is the dual space of H1(Ω). THEOREM 2.1. The continuous problem (2.6) with ǫ= 0 satisfies the following energy equality: ku(t)k2+λk∇d(t)k2+ 2 Zt 0 (νk∇u(s)k2+λ γk∂td(s) + (u(s)· ∇)d(s)k2)ds =ku0k2+λk∇d0k2+ 2 Zt 0 hg(s),u(s)ids,(2.7) that holds for any t∈[0, T ]. On the other hand, for ǫ > 0, the system satisfies ku(t)k2+λk∇d(t)k2+ 2 Zt 0 (νk∇u(s)k2+λ γk∂td(s) + (u(s)· ∇)d(s)k2)ds +λǫ2 2kq(t)k2=ku0k2+λk∇d0k2+λǫ2 2kq0k2+ 2 Zt 0 hg(s),u(s)ids,(2.8) for any t∈[0, T ], where q0= (|d0|2−1)/ǫ2. Proof. Defining w(d,u) := ∂td+ (u· ∇)dthe following identity holds: λ γc(u,d,w) = λ γkwk2−λ γ(∂td,w).(2.9) On the other hand, we can re-write Eq. (2.6a) as: (w,¯ d) + γ(∇d,∇¯ d) + γbd(q, d,¯ d) = 0. Taking ¯ d=∂td, we easily get: (2.10) −(w, ∂td) = γ 2∂tk∇dk2+γ bd(q, d, ∂td). The time derivative of Eq. (2.6b) leads to 2bd(¯q, ∂td,d) + ǫ2(∂tq, ¯q) = 0, which for ¯q=qallows to write the last term of (2.10) as follows bd(q, d, ∂td) = ǫ2 2(q, ∂tq) = ǫ2 4∂tkqk2. 6 Hence, we finally have that (2.9) is expressed as: λ γc(u,d,w) = λ γkwk2+λ 2∂tk∇dk2+λǫ2 4∂tkqk2. The desired energy equality is obtained testing (2.6c)-(2.6d)against (u, p), using the previous equality and the skew-symmetry property c(u,¯ u,¯ u) = 0 for any ¯ u∈H1 0(Ω). We get: ∂tkuk2+λk∇dk2+λǫ2 2kqk2+ 2 νk∇uk2+λ γkwk2= 2 hg,ui a.e. in (0, T ). For ǫ= 0, the energy equality (2.7) is obtained after integrating the previous equation in the time interval (0, t). The energy equality for the penalized problem (2.8) is proved by noting that Eq. (2.6b) holds at t= 0 with q0= (|d0|2−1)/ǫ2for ǫ > 0. Pressure stability relies on the well-known inf-sup condition: inf ¯p∈L2 0(Ω) sup ¯ u∈H1(Ω) bu(¯p, ∇ · ¯ u) k¯pkk¯ uk1 ≥βu>0, which is known to be true due to the surjectivity of the divergence operator from H1 0(Ω) onto L2 0(Ω);L2 0(Ω) is the space of L2(Ω) functions with zero mean value. We refer to [31, 18] for some regularity results for the pressure in the transient Navier-Stokes system. The control over the Lagrange multiplier qis not so well understood. The GL penalty version introduces L2(Ω) control over qthat is lost when the penalty ǫց0. So, this stability does not apply for the singular limit ǫ= 0, and the well-posedness of the original problem can only rely on an inf-sup condition. The inf-sup condition inf ¯q∈H−1(Ω) sup ¯ d∈H1 0(Ω) bd(¯q, d,¯ d) k¯qk−1k¯ dk1 ≥βd(d)>0,(2.11) has been proved very recently in [19] under some regularity assumptions over din the frame of the steady harmonic maps problem, which is (2.6a)-(2.6b) without the time derivative and convective terms. This regularity is much stronger than the one that the energy estimate (2.7) provides for d. With regard to the long-term behavior for a zero forcing term (g= 0), we can see that for t→ ∞ the energy dissipation in (2.7) goes to zero: νk∇u(t)k2+λ γk∂td(t) + (u(t)· ∇)d(t)k2→0, hence ugoes to the trivial stationary point u= 0 and k∂tdk → 0. The dcomponent of the solution exhibits non-trivial stationary states. Such a stationary states are minima of the Oseen-Frank free-energy function (1.1). It means that there exist steady solutions with k∇dk>0. This fact is shared by both the limit and penalized case. For the penalized case, there is an additional term in the energy, which is the penalty energy λǫ2 2kqk2. 7 3. Spatial discretization. Let Thbe a partition of Ωinto a set of finite elements {K}. For every element K, we denote by hKits diameter, and set the characteristic mesh size as h= maxK∈ThhK. The space of polynomials of degree less or equal to k > 0in a finite element Kis denoted by Pk(K). The space of continuous piecewise polynomials is defined as (3.1) Pk h=vh∈ C0(Ω) such that vh|K∈ Pk(K)∀K∈ Th. These approximations are usually called H1-conforming approximations, because of the inter-element continuity. The space Pk his spanned by the set of nodal functions {πa h}a∈Nh, where Nhis the set of nodes in the mesh. Therefore, any function ϕh∈ Pk hcan be uniquely determined in terms of its nodal values {ϕa}a∈Nhas Pa∈Nhϕaπa h(see [5, 13]). The nodal interpolation of a continuous function ϕ∈ C0(Ω) is denoted by πh(ϕ) = Pa∈Nhϕ(xa)πa h Let us consider a conforming finite element discretization of problem (2.6). The finite element space Dhfor the director vector dhis chosen to be (P1 h)d. We also consider the space Qhfor the Lagrange multiplier qhto be P1 h. The constraint form reads as bd(qh,dh,¯ dh) = (qh,dh·¯ dh).(3.2) In this case, the constraint is satisfied in a discrete sense, as the incompressibility condition for the fluid problem. In the frame of harmonicmaps, Hu et al. have considered the following modification of the constraint form in [19]: bd(qh,dh,¯ dh)h= (qh, πh(dh·¯ dh)).(3.3) Since πh(·)∈Qhwith the previous choice of finite element spaces, the constraint equation amounts to saying that |da|= 1 for any a∈ Nh. Furthermore, the finite element pair Dh×Qhhas been recently proved to satisfy the corresponding discrete version of the inf-sup condition (2.11): inf qh∈Qh sup ¯ dh∈Dh bd(qh,dh,¯ dh)h kqhk−1k¯ dhk1 ≥βh,d(dh)>0,(3.4) in [19]. It has allowed to provethe well-posedness of the tangent problem that arises from the full Newton linearization of the steady-state harmonicmaps problem, in the vicinity of a local minimum under strong regularity assumptions. In this proof, the fact that the projection πh(·) has been used is necessary, so it does not apply to bd(qh,dh,¯ dh). In any case, numerical experimentation says that this choice is also stable for the problems considered in Section 6. REMARK 3.1. Let us recall that, for Dirichlet boundary conditions over dh, the discrete inf-sup condition (3.4) is violated when there are elements K∈ Thwith all the nodes constrained; we can easily check that there is no control of the Lagrangemultiplier at a boundary node that is connected to boundary nodes only. This type of meshes usually leads to problems and should be avoided. In any case, an alternative way to circumvent this problem is to 8 consider homogeneous Dirichlet boundary conditions over qtoo, as proposed in [19]. Since |da|= 1 for appropriate boundary data, qa= 0 is an appropriate condition. In [19], the saddle-point version of the harmonic maps problem is used together with a Newton linearization, and the corresponding inf-sup condition for the tangent problem is proved. As an alternative, it could be used for the minimization step of the linearized problem in the popular Alouges’ method proposed in [1] (see also [3, 2] and Section 5.3). In the frame of liquid crystals, we consider two different solvers, that are extensions of these two approaches to the problem at hand. Let us point out that the dh−qhblock matrices in the corresponding linear system are diagonal matrices when using closed (nodal)integration, and so computationally more efficient; closed (nodal) integration of the constraint terms is the way to get a lumped mass matrix for Lagrangian elements. The closed integration of the constraint trilinear form consists of: bd(qh,dh,¯ dh)hs =X a∈Nh qada·¯ daZΩ πa hdx= (qh,dh·¯ dh)s,(3.5) where the sub-label sindicates lumped sub-integration of the term. The nodal enforcement of the restriction is even more explicit in this case. Here and in the sequel, we will also use the sub-index sfor inner-products (f, g)s=RΩπh(fg) dxand L2-discrete norm kfks= (f, f)1/2 sinvolving finite element functions to denote that closed integration is used. REMARK 3.2. The nodal enforcement of the constraint could also be understood as a collocation method for the constraint equation. In this case, the discrete version of the Lagrange multiplier space consists of Qh={δ(xa), a ∈ N}, where δ(xa) : C0(Ω) →R is defined by hδ(xa), vi=v(xa)for v∈ C0(Ω). This approach to the problem is not so powerful, because it can only be used for the limit case, since the penalty term is ill-posed. For the Navier-Stokessub-problem we consider the standard MINI element, in which the pressure finite element space Phis taken as P1 h, and the velocity space Vhis (P1 h)d⊕(Bh)d, where Bh={vbsuch that vb|K∈ P3(K), vb|∂K = 0, vb|K≥0,∀K∈ Th} is the space of bubbles (cubic in dimension 2) at every element (see e.g. [7]). This velocitypressure finite element pair is known to satisfy the discrete inf-sup condition inf ¯ph∈Qh sup ¯ uh∈Vh bu(¯ph,¯ uh) k¯phkk¯ uhk1 ≥βu,h ≥0 for βu,h uniform with respect to h. Onwards, we consider the skew-symmetric form (for Dirichlet boundary conditions) of the convective term in the uhproblem: ec(uh,vh,¯ uh) = ((uh· ∇)vh,¯ uh) + 1 2((∇ · uh)vh,¯ uh). 9 The finite element approximation of system (2.6) reads as: find (dh(t), qh(t),uh(t), ph(t)) ∈ Dh×Qh×Vh×Phsuch that (∂tdh,¯ dh) + γ(∇dh,∇¯ dh) + c(uh,dh,¯ dh) + γbd(qh,dh,¯ dh)δρ = 0,(3.6a) bd(¯qh,dh,dh)δρ −ǫ2(qh,¯qh)ρ= (1,¯qh)ρ,(3.6b) (∂tuh,¯ uh) + ν(∇uh,∇¯ uh) + ec(uh,uh,¯ uh) +λ γc(¯ uh,dh, ∂tdh+ (uh· ∇)dh)−bu(ph,¯ uh) = hg,¯ uhi,(3.6c) bu(¯ph,uh) = 0,(3.6d) at almost every t∈(0, T ], for any (¯ dh,¯qh,¯ uh,¯ph)∈Dh×Qh×Vh×Ph. The sub-label δ takes the values hwhen using πh(·)in bd;ρtakes the value swhen sub-integration is used. REMARK 3.3. Many existing liquid crystal finite element approximations involve an auxiliary variable (see e.g. [17, 4]). We can re-formulate the problem, by introducing an auxiliary variable whand its corresponding finite element space Wh; Eq. (3.6a) is replaced by (wh,¯ dh) + γ(∇dh,∇¯ dh) + γ bd(qh,dh,¯ dh)δρ = 0, (∂tdh,¯ wh) + c(uh,dh,¯ wh)−(wh,¯ wh) = 0, for any (¯ dh,¯ wh)∈Dh×Wh. Then, the elastic stress term c(¯ uh,dh, ∂tdh+ (uh· ∇)dh)is replaced by ZΩ (∇dh)twh·¯ uhdx. This approach would in principle introduce extra unknowns to the problem, which would preferably be avoided. In most existing liquid crystal algorithms, always for the penalized GL problem, the space for Whis taken equal to Dh, i.e. (P1 h)d(see [17, 4]). The method is proved to be stable and convergent but, since (uh·∇)dh/∈Dh, the problem for whis global and whcannot be locally eliminated. However, taking Wh≡(Vh· ∇)Dhand Dh⊂Wh, we simply have wh=∂tdh+ (uh· ∇)dh. Implicitly, this is the approach that has been used above and the one used in [23] for the approximation of the penalized GL method. The saddle-point approach has a clear advantage with respect to the previous finite element approximations of liquid crystal flows. In the frame proposed herein, we can approximate numerically both the limit problemand GL penalized versions using the same numerical approximation. Existing algorithms [23, 25, 26, 4, 17] could not take ǫarbitrarily small, since the condition number of the matrix blows up with ǫ−2. As a rule of thumb, ǫ= 0.05 was the limit value used in numerical experiments. Since the stability of the saddle-point structure as ǫց0is kept by virtue of the discrete inf-sup condition (3.4), the linear system is non-singular even for the limit problem. In the numerical experiments section we analyze 16 5.1. Exact Newton scheme. It is clear that both operators are nonlinear. At this point, we can linearize the problem using an exact Newton linearization. Given a previous iterate (dn+1,k h, qn+1,k h,un+1,k h, pn+1,k h), the new iterate (dn+1,k+1 h, qn+1,k+1 h,un+1,k+1 h, pn+1,k+1 h) = (dn+1,k h, qn+1,k h,un+1,k h, pn+1,k h) + (δdk+1 h, δqk+1 h, δuk+1 h, δpk+1 h) is obtained after solving the linear system: *dLd(dn+1,k h, qn+1,k h,un+1,k h) d(dh, qh,uh)·(δdk+1 h, δqk+1 h, δuk+1 h),(¯ dh,¯qh)+ =−DLd(dn+1,k h, qn+1,k h,un+1,k h),(¯ dh,¯qh)E,(5.2) *dLu(un+1,k h, pn+1,k h,dn+1,k h) d(uh, ph,dh)·(δuk+1 h, δpk+1 h, δdk+1 h),(¯ uh,¯ph)+ =−DLu(un+1,k h, pn+1,k h,dn+1,k h),(¯ uh,¯ph)E, where dF(x∗) dx·δx∈Y′denotes the weak Gˆateaux derivative of the functional F:X→Y at x∗with respect to xin the direction δx∈X, for X,YBanach spaces. Problem (5.2) is the tangent problem of (5.1) around (dn+1,k h, qn+1,k h,un+1,k h, pn+1,k h). Every nonlinear iteration kof the exact Newton linearization requires to solve a linear system coupling all the unknowns (δdk+1 h, δqk+1 h, δuk+1 h, δpk+1 h)in the problem, with the corresponding computational cost; in dimension 3, it involves eight degrees of freedom per node. 5.2. Quasi-Newton scheme. For numerical purposes, it is convenient to decouple the differentvariables in the problem, in orderto reduce the CPU time and memory usage. In particular, since we want to decouple sub-problems (4.2a)-(4.2b) and (4.2c)-(4.2d), we consider a quasi-Newton method in which the tangent matrix decouples problems. Given the iterate (dn+1,k h, qn+1,k h,un+1,k h, pn+1,k h), the new iterate is obtained after solving the two linear sub-problems: *dLd((dn+1,k h, qn+1,k h),un+1,k h) d(dh, qh)·(δdk+1 h, δqk+1 h),(¯ dh,¯qh)+ =−DLd((dn+1,k h, qn+1,k h),un+1,k h),(¯ dh,¯qh)E, *dLu((un+1,k h, pn+1,k h),dn+1,k+1 h) d(uh, ph)·(δuk+1 h, δpk+1 h),(¯ uh,¯ph)+ =−DLu((un+1,k h, pn+1,k h),dn+1,k+1 h),(¯ uh,¯ph)E, Doing that, we have reduced the number of terms in the tangent problem, since we have neglected cross derivatives, i.e. the variation of Ldwith respect to uhand vice-versa. A nice property of this approach is the fact that it allows modularity. Separate codes can be 17 used to solve the two sub-problems, only needing to pass the unknowns from one code to the other at every iteration. So, even though this method involves nonlinear iterations, it has been linearized in such a way that the computation of (uh, ph)is segregated from the one for (dh, qh), notably reducing the solver CPU time per iteration with respect to the semi-implicit method. 5.3. Nonlinear block Gauss-Seidel scheme with Alouges’ method. The linearization of equations (4.2a)-(4.2b), that are equivalent to the harmonic maps problem plus a convection term, can be linearized by extending the strategy proposed by Alouges in [1]. We refer to [1, 3] for a mathematical analysis of the Alouges’ method in the frame of steady harmonic maps. The idea is to consider a Picard linearization of (4.2a)-(4.2b) together with a normalization of the director field approximation. Obviously, this approach only has sense when the constraint over dhis exactly enforced on the nodes, using one of the alternatives proposed above, and no penalty is introduced, i.e. ǫ= 0. Let us consider the previous iterate (dn+1,k h, qn+1,k h,un+1,k h, pn+1,k h). First, we compute (˜ dn+1,k+1 h, qn+1,k+1 h), solution of the linear system (δt˜ dn+1,k+1 h,¯ dh) + γ(∇˜ dn+1,k+1 h,∇¯ dh) +c(un+1,k h,˜ dn+1,k+1 h,¯ dh) + γbd(qn+1,k+1 h,dn+1,k h,¯ dh)hρ = 0, bd(¯qh,dn+1,k h,˜ dn+1,k+1 h)hρ =h1,¯qhiρ. Then, we compute dn+1,k+1 h, as the normalization of ˜ dn+1,k+1 hon the nodes. So, at every node, we compute (dn+1,k+1)a=(˜ dn+1,k+1)a |(˜ dn+1,k+1)a|. Obviously, dn+1,k+1 hsatisfies the nonlinear constraint bd(¯qh,dn+1,k+1 h,˜ dn+1,k+1 h)hρ =h1,¯qhiρ. The third step of the algorithm consists of solving (4.2c)-(4.2d) with the known value of dn+1,k+1 h, e.g. using a Picard linearization: we compute (un+1,k+1 h, pn+1,k+1 h)solution of (δtun+1,k+1 h,¯ uh) + ν(∇un+1,k+1 h,∇¯ uh) + ec(un+1,k h,un+1,k+1 h,¯ uh)h+bu(pn+1,k+1 h,¯ uh) +λ γc(¯ uh,dn+1,k+1 h, ∂tdn+1,k+1 h+ (un+1,k+1 h· ∇)dn+1,k+1 h) = gn+1,¯ uh, bu(¯ph,un+1,k+1 h) = 0. So, the final procedure involves a linearized harmonic map-like system, with a convection term, a normalization of the director field, and the linearized Navier-Stokes equations, e.g. using Picard, evaluating the coupling elastic term with the director field of the second step. The final problem has a computational cost per iteration similar to the quasi-Newton algorithm above, since the two sub-problems have been decoupled. 18 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 60 70 80 Time Energy Implicit Implicit P1 Semi−implicit S−imp small dt (a) Ed(t) 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 1.5 2 Time Energy Implicit Implicit P1 Semi−implicit S−imp small dt (b) Eu(t) FIG. 6.1. Example 1: Ed(t)and Eu(t)plots for the nodal implicit, P1implicit and explicit methods, for ǫ= 0 and k= 10−2. The explicit method is also used with k= 10−3. 0 0.002 0.004 0.006 0.008 0.01 5 5.5 6 6.5 7 7.5 8 log(cond. number) ε h = 2.0e2 h = 2.0e3 h = 2.0e4 h = 2.0e5 (a) Cond. number vs. ǫ −2 −1.5 −1 −0.5 4 5 6 7 8 9 10 log(cond. number) log(h) ε = 0.0 ε = 1e−2 ε = 1e−3 slope 4 (b) Cond. number vs. h FIG. 6.2. Example 1: Condition number of the dblock-matrix vs. (ǫ, h). −2 −1.5 −1 −0.5 2 3 4 5 6 7 log(norm) log(h) H−1 norm L2 norm H1 norm FIG. 6.3. Example 1: qhstability vs. hfor the L2,H1and discrete H−1norm. 6. Numerical experimentation. In this section, we perform some numerical experiments, in order to check the behavior of the methods proposed above. We will distinguish between three different numerical methods: •The implicit method with nodally exact enforcement of the constraint and Crank- 19 Nicolson time integration. Thus, the method consists of system (4.2) with ρ=s and δ=h, that is to say, using a closed integration rule for the constraint equation and the bilinear form bdin (3.3), i.e. Eq. (3.5). The problem is linearized using the quasi-Newton scheme in Section 5.2, that decouples dhand uhcomputation at the linear solver level. We will denote this method as nodal implicit method. •The implicit method with P1 has Lagrange multiplier space and Crank-Nicholson time integration. In this case, the method consists ofsystem (4.2) with the expression of bdin (3.2). Again, we use the quasi-Newton scheme for linearization. We will denote this method as P1implicit method. •The semi-implicitmethod (4.9) with a closed integration rule for the constraint equation and the bilinear form bdin (3.3). We will denote this method as semi-implicit method. For all the methods, we will consider both exact and penalized formulations. In all cases, we have used quadrature rules that integrate exactly all the terms in the linear system. One of the outputs of the simulations are the time behavior of the different energies interacting in the system. Let us define the elastic and kynetic energies respectively as: Ed(t) = k∇dh(t)k2,Eu(t) = kuh(t)k2. The validationof the code has been carried out by using the method of manufacturedsolutions and checking that the implicit methods under consideration are energy preserving, i.e. the equality (4.4) is satisfied up to convergence tolerance. So, for the implicit methods with ǫ= 0 and zero forcing terms, the energy equality has been checked: λEd(tn) + Eu(tn) + n−1 X m=0 Edis(tm+1 2) = E0 where Edis(t) = 2νkk∇uh(t)k2+ 2λ γkkδtdh(t) + uh(t)· ∇dh(t)k2 denotes the energy dissipation and E0the initial energy (see (4.4)). We have analyzed three different test problems. Two of them have been presented in previous numerical works about the approximation of liquid crystals and serve for comparison with pre-existing techniques. The third example is a problem with known analytical solution that clearly serves to assess the algorithms. 6.1. Example 1: A smooth harmonic map. The first test under consideration has been previously solved in [4]. We consider the problem (2.2) for ǫ= 0 in the square domain Ω = (−1,1)2and g=0. The initial conditions are: u0=0,d0= (sin(a),cos(a))t, a = 2.0π( cos(x)−sin(y)), 20 whereas homogeneous Dirichlet and Neumann boundary conditions are enforced over uand drespectively. So, the only energy introduced to the system is via the initial energy E0. The physical parameters are λ=γ= 1.0and ν= 0.1, and the numerical parameters are chosen as h= 2−5,k= 0.01, and ǫ= 0, unless otherwise stated. In Fig. 1 we show the elastic and kynetic energies Ed(t)and Eu(t)with respect to time. The initial conditions introduce energy to the system via the initial elastic energy λk∇d0k2, that is related to (4.5) (see [4, Fig. 5.1]). It is clear out of the results that most of the initial elastic energyis (unsurprisingly)transferredto Ed, whereas the kynetic energyfor the velocity Euat its peak (t= 0.1 s) is only about a 5% of Ed. The steady-state that is reached for the initial condition stated above has no energy as t→ ∞, since the steady-state dhis constant in space (see [4, Fig. 5.1]). We have plot the energy results for the three methods under consideration, whereas for the semi-implicit method, we have also included the results for k= 10−3. Out of these plots, we can conclude that the nodal constraint and the P1Lagrange multiplier lead to very similar results. As expected, the semi-implicit method shows some lag in the dynamics to the steady state, but reducing by 10 the time step size, the results are almost identical to those with an implicit method. Now, we consider λ= 0.0, reducing the problem to the transient harmonic maps system. We do this since we want to evaluate the dependence of the condition number and the constant in the discrete inf-sup condition (3.4), for the dproblem; these features are well-known for the Navier-Stokes block. We compute the condition number of the system matrix for a given time step value, a penalty parameter ǫ= 0,10−2,10−3and h= 2−iwith i= 2,3,4,5; the results are collected in Fig. 6.2. In 6.2(a) we see that the condition number of the matrix is not sensitive to ǫ, which means that the stability of the Lagrange multiplier comes from the discrete inf-sup condition. Let us remark that for the GL formulations in the literature (see e.g. [4, 23, 17, 25, 26]) the condition number blows up for ǫց0. We also plot in Fig. 6.2(b) the condition number in terms of the mesh size hfor different values of ǫ. We see that the condition number is almost O(h−4). Finally, we want to analyze the constant in the inf-sup condition (3.4). In order to do this, we evaluate the L2norm, the H1norm and a discrete H−1norm. Given qh∈Qh, let us find the discrete Riestz projection gh∈Qhsuch that (gh, fh)H1=hqh, fhiH−1×H1. We define the discrete H−1norm as follows: kqhkH−1 d:= sup fh∈Qh hqh, fhiH−1×H1 kfhk1 = sup fh∈Lh (gh, fh)H1 kfhk1 =kghk1. In Fig. 6.3 we plot the norm of the Lagrange multiplier at a given time step size, for the uniform partition introduced above. We can easily see that the Lagrange multiplier qhis stable in both the H−1and L2norms. As expected, the H1norm blows up with h. This test indicates that for smooth enough solutions, qh∈L2(Ω). 21 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 Time Energy Semi−implicit Nodal implicit P1 implicit (a) ǫ= 0.100 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 Time Energy Semi−implicit Nodal implicit P1 implicit (b) ǫ= 0.075 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 Time Energy Semi−implicit Nodal implicit P1 implicit (c) ǫ= 0.050 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 Time Energy Semi−implicit Nodal implicit P1 implicit (d) ǫ= 0.025 FIG. 6.4. Example 2: Elastic energy Ed(t)plots for the nodal implicit, P1implicit and semi-implicit methods. The results are shown for different values of ǫ. 6.2. Example 2: Annihilation and stable defects. The second test we consider can be found in any article about GL-based numerical approximations of liquid crystals. It consists on the annihilation of regularized initial singularities. The problem is solved in the square domain Ω = (−1,1)2and g=0. The initial conditions are: u0=0,d0=˜ d/q|˜ d|2+ǫ2,for ˜ d= (x2+y2−0.25, y)t. The physical parameters are λ=γ=ν= 1.0, and the numerical parameters are chosen as h= 2−5,k= 0.001, unless otherwise stated. We have considered ǫ= 0,2.5·10−3,5· 10−3,7.5·10−3,10−2. As commented in [4, Section 5.2] the solution critically depends on the initial data and so, on ǫ. It was also pointed out by Lin and Liu in [23, Example 2.2] that the evolution of the singularities was very sensitive to the algorithm and the mesh selected in the computation. Furthermore, in both references the authors point out that the local energy at the singularity goes to infinity as ǫց0. In this article, we put particular attention on the solution with respect to ǫ. The initial condition introduces initial energy to the system in two different ways. On one side, via the initial elastic energy, and on the other side, via the q0term. 22 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 Time Energy Semi−implicit Nodal implicit P1 implicit (a) ǫ= 0.100 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 Time Energy Semi−implicit Nodal implicit P1 implicit (b) ǫ= 0.075 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 Time Energy Semi−implicit Nodal implicit P1 implicit (c) ǫ= 0.050 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 Time Energy Semi−implicit Nodal implicit P1 implicit (d) ǫ= 0.025 FIG. 6.5. Example 2: Kynetic energy Eu(t)plots for the nodal implicit, P1implicit and semi-implicit methods. The results are shown for different values of ǫ. However, the second term reads as: λǫ2 2kq0k2=λ 2kd0·d0−1k2=λ 2ZΩǫ2 |˜ d|2+ǫ22 ,(6.1) which tends to zero as ǫց0. At the discrete level the energy also depends on the projection that is used for obtaining d0 hfrom d0. We use a nodal projection for the sake of comparison with previous works. However, this projector is ill-posed as (ǫ, h)ց0in the singularity, but can be solved by giving an average value for the node over the singularity, if there is one. In the limit case (ǫ= 0), the initial energy is only coming from the initial elastic energy. In Fig. 6.4 we show the elastic energy Ed(t)for the three formulations and different values of ǫ. The implicit saddle-point Ginzburg-Landau formulation proposed herein leads to very similar results to those obtained for classical GL formulation (see [17, 4, 23]). On the other hand, it is clear that the annihilation time goes to ∞as ǫց0. For ǫ= 0.025 the implicit methods converge to a stable solution with non-zero elastic energy and zero velocity. In fact, as commented above, this long-term behavior is in agreement with the structure of the dynamical system at hand. The steady-state solution for dhis a solution of the discrete harmonic maps problem. These numerical results are in concordance with the mathematical 23 (a) t= 0.00 (b) t= 0.25 (c) t= 0.30 (d) t= 1.00 FIG. 6.6. Example 2: Director vector dfields at different time step values for the nodal implicit method and ǫ= 5 ·10−3. analysis of Brezis, Coron and Lieb in [6], which proves that this singular solution minimizes the energy of a stationary pair of point defects. A physical justification of the existence of stable combed hedgehog point defects can be found in [28] using different approaches, e.g. Peach-K¨ohler type arguments (see also [30, 29, 14, 9]). The results obtained with nodal and P1implicit methods are almost identical. We show the director vector fields at different time values for ǫ= 5 ·10−3in Fig. 6.6. We see how the singularities approach and annihilate at some point over t= 0.30 s, which coincides with an abrupt decreasement of elastic energy (see Fig. 6.4(c)). In Fig. 6.4 we also show the results for the explicit method. It is clear, out of these results, that the effective penalty of the semi-implicit formulation is smaller that the one for the implicit method. For ǫlarge enough or small enough, the results are fairly similar. However, in the transition between them, the explicit method exhibits a lag with respect to the implicit solution, or can even tend to a different stationary point. Summarizing, the solution is critically ǫ-dependent. For ǫ= 0.0, the results are almost identical to those for ǫ= 0.025. As far as we know, this behavior for ǫց0has not been previously analyzed because previous GL formulations were ill-posed in this assymptotic regime, and a minimum value of ǫ= 0.05 was common practice. Herein, this problem has been solved by using a saddle-point GL approach. 24 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 Time Energy λ = 1.0 λ = 0.0 (a) Ed(t)vs. λ 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 Time Energy λ = 1.0 λ = 0.0 (b) Eu(t)vs. λ 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 Time Energy h = 1/32 h = 1/48 (c) Ed(t)vs. h 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 Time Energy h = 1/32 h = 1/48 (d) Eu(t)vs. h 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 Time Energy k = 0.010 k = 0.001 (e) Ed(t)vs. k 0 0.2 0.4 0.6 0.8 1 0 0.01 0.02 0.03 0.04 0.05 0.06 Time Energy k = 0.010 k = 0.001 (f) Eu(t)vs. k FIG. 6.7. Example 2: Ed(t)and Eu(t)plots for the nodal implicit method with ǫ= 5 ·10−3and different values of (λ, h, k). We also plot the kynetic energy in all cases (see Fig. 6.5). These results make clear that this test problem is almost a harmonic maps problem, since the energy transferred to uis very small. The results for the implicit methods are very similar, but for ǫ= 2.5·10−3the implicit P1method exhibits two picks of kynetic energy whereas the nodal version only one. In Fig. 6.7(a)-6.7(b)we have considered the full liquid crystal problem with λ= 1.0and the transient harmonic maps problem, i.e. λ= 0.0. The elastic energy is almost identical in 25 0 0.2 0.4 0.6 0.8 1 20 25 30 35 40 45 50 Time Energy w = 50.0 w = 5.0 FIG. 6.8. Example 2: Es(t)plot for the nodal implicit method, ǫ= 0.0and ω= 5,50. (a) t= 0.04 (b) t= 0.08 (c) t= 0.12 (d) t= 0.23 FIG. 6.9. Example 2: Director vector dfields at different time step values for the nodal implicit method, ǫ= 0 and ω= 50. both cases, and so, the uinfluence almost neglectable. Obviously, the kynetic energy is zero for λ= 0.0in Fig. 6.7(b). With regard to the mesh size, we have compared both elastic and kynetic energy for h= 2−5,2−6, showing that the solution is well converged for the mesh being used. Analogously, we show the results for k= 10−3,10−2, that allows to say that the time step size is also acceptable. Since the initial energy introduced to the system is not enough for the annihilation of