An overview on numerical analyses of nematic liquid crystal flows

An overview on numerical analyses of nematic liquid crystal flows S. Badia∗F. Guill´en-G´onzalez †J. V. Guti´errez-Santacreu ‡ June 19, 2010 Abstract The purpose of this work is to provide an overview of the most recent numerical developments in the field of nematic liquid crystals. The Ericksen-Leslie equations govern the motion of a nematic liquid crystal. This system, in its simplest form, consists of the Navier- Stokes equations coupled with an extra anisotropic stress tensor, which represents the effect of the nematic liquid crystal on the fluid, and a convective harmonic map equation. The sphere constraint must be enforced almost everywhere in order to obtain an energy estimate. Since an almost everywhere satisfaction of this restriction is not appropriate at a numerical level, two alternative approaches have been introduced: a penalty method and a saddle-point method. These approaches are suitable for their numerical approximation by finite elements, since a discrete version of the restriction is enough to prove the desired energy estimate. The Ginzburg-Landau penalty function is usually used to enforce the sphere constraint. Finite element methods of mixed type will play an important role when designing numerical approximations for the penalty method in order to preserve the intrinsic energy estimate. The inf-sup condition that makes the saddle-point method well-posed is not clear yet. The only inf-sup condition for the Lagrange multiplier is obtained in the dual space of H1(Ω). But such an inf-sup condition requires more regularity for the director vector than the one provided by the energy estimate. Herein, we will present an alternative inf-sup condition whose proof for its discrete counterpart with finite elements is still open. 2000 Mathematics Subject Classification. 35Q35, 65M12, 65M60. Keywords: liquid crystal, Navier-Stokes, stability, convergence, finite elements, penalization. ∗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. E.D.A.N., University of Sevilla, Aptdo. 1160, 41080 Sevilla, Spain. E-mail: [email protected]. This author’s work was partially supported by project MTM2009-12927, Spain. ‡Dpto. de Matem´atica Aplicada I, University of Sevilla, E. T. S. I. Inform´atica. Avda. Reina Mercedes, s/n. 41012 Sevilla, Spain. [email protected]. This author’s work was partially supported by project MTM2009-12927, Spain. 1 1 Introduction Liquid crystals are commonly considered as the fourth state of matter, different to gases, liquids, and solids. This is due to the fact that liquid crystals exhibit phases between a liquid and a crystalline solid which are known as mesophases. There are various types of liquid crystals, according to the degree of positional or orientational ordering shown by the molecules that compose them. Different degrees of ordering can be achieved, depending on the temperature (thermotropic) and/or the concentration of a solute in a solvent (lyotropic). The simplest liquid crystal phase is the nematic one, which is made of elongated rod-like molecules with similar size, whose centers of mass have no positional order (as in an isotropic liquid), but tending to align along certain locally preferred directions, confering the anisotropic structure. The orientational order is typically modeled by a unit vector field, the director field d, with |d|= 1, which represents the average orientation of the long axes of the molecules in a volume element at a point. As the temperature is lowered in the thermotropic case, or the concentration is increased in the lyotropic case, new thermodynamical states appear, and the molecules begin to separate into parallel, equally spaced layers. This is the smectic liquid crystal phase. In a smectic liquid crystal, the positions of the molecular centers of mass flow freely in each layer without correlation from one layer to the next. The angle between the director field dand the local smectic layer normal nis denoted by θ. When θ= 0, the structure is called a smectic-A liquid crystal, otherwise it is called a smectic-C liquid crystal. On some occasions there are transitions from smectic-A to the smectic-C phase, i.e. the angle θgrows smoothly from zero in response to a decrease of temperature or an increase of concentration. In the last years, the study of liquid crystals has aroused an increasing interest in biology (cell membranes), physics and engineering (in the growing technological industry of electronic devices) owing to their optical properties. The alignment of the director field dis affected by applied electric or magnetic fields, that can rotate the director so that it is aligned parallel to them; the molecules of a liquid crystal exhibit dielectric or diamagnetic properties. The director field can also be affected by boundary conditions. A usual geometry for liquid crystal display (LCD) devices is that of a thin film made of pixels which are filled by a liquid crystal. Each pixel consists of two transparent electrodes, and two polarizing filters. The surfaces are treated to keep the director of molecules at the top and the bottom perpendicular to each other and to the surface normals. The molecules in the pixel describe a helicoidal structure in the transition between surfaces. If a polarized ray of light parallel to the molecules in the external surface enters the pixel, the plane of the polarization rotates with the director field. Under external electric fields normal to the surfaces, the helicoidal structure is broken, the rotation is not possible, and light is blocked by the polarizer. On the 2 other hand, liquid crystals also provide a description of some interesting materials, as DNA or petroleum. The two main phenomenological theories describing spatial configurations in nematic liquid crystals are the Oseen-Frank [53, 25] and Landau-de Gennes [20, 21] theories. Both approaches consist in modeling equilibrium states as minima of a free-energy functional. Such functionals are constructed subject to symmetry and invariance principles, to capture some properties observed from experiments. The Oseen-Frank free energy is considered as a functional of the director vector d. In its most basic form, the free energy functional is given by E(d) = ZΩ {K1|∇ · d|2+K2(d·(∇ × d))2+K3(d×(∇ × d))2},(1) where K1,K2, and K3are the splay,twist, and bend elastic constants, respectively. Note that when these constants are equal, the Dirichlet energy becomes E(d) = KZΩ |∇d|2. Upon minimizing this energy subject to the sphere constraint |d|= 1, the following optimality system appears −∆d− |∇d|2d= 0 in Ω.(2) The Oseen-Frank theory is limited in the sense that it can only explain point defects in liquid crystal materials but not the more complicated line and surface defects that are also observed experimentally. The defect points or singularities in liquid crystals are regions where the anisotropic properties of molecules are broken. That is, the liquid crystal behaves as an isotropic fluid. Therefore, the director field cannot be defined. Mathematically, they are modeled by |d|= 0. The defect points can be achieved by means of the boundaries conditions. The Landau-de Gennes functional is established in terms of the order parameter tensor Q (traceless and symmetric) as E(d) = ZΩ {K1 2|∇Q|2+K2 2|∇ · Q|2+atr(Q2)−btr(Q3) + c(tr(Q2))2}, where Qdepends on the components of the director vector das Qij =didj−1 dδi,j;δis the Kronecker delta, dis the space dimension and tr(·) indicates the trace of the matrix. This theory was one of the major achievements of P. G. de Gennes, who was awarded a Nobel prize in physics in 1991. An alternative strategy to study the motion of defect points in liquid crystals is to consider the long-time behavior of the harmonic map flow for which it is also interesting to incorporate the influence of the velocity. On the contrary, in many situations, the anisotropic local orientation of 3 the director field influences the stress tensors that govern the fluid velocity. The hydrodynamic theory of nematic liquid crystals was established by Ericksen [22, 23] and Leslie [39, 40]. The fundamental system consists of a set of fully coupled, macroscopic equations, that contains the Oseen-Frank elastic theory governing the steady state, equilibrium solutions. A variant of the Ericksen-Leslie equations was proposed by Lin in [41]. The equation therein is written in such a way that the sphere constraint is not explicitly enforced. The obtention of the energy law for this equation uses the fact that the director field satisfies the constraint everywhere. So, it is not suitable for numerical purposes, since the numerical approximation does not satisfy the constraint everywhere, and so, does not have an associated energy law. Then, some alternatives are introduced. In order to get numerical methods with an associated energy law, the sphere constraint is usually penalized with the Ginzburg-Landau penalty function such as in the works of Becker, Feng and Prohl [8], and Lin, Liu, and Zhang [45]. An alternative recently proposed by Badia, Guill´en-Gonz´alez and Guti´errez-Santacreu in [5] is to use an equivalent saddle-point formulation of the system proposed by Lin in [41]. It provides a system of partial differential equations equivalent to the one in [41] that also leads to numerical methods with an energy law. Furthermore, a Ginzburg-Landau-type penalization can be introduced in the saddle-point version, and treat both the original and penalized problems using a unified numerical approach. The goal of this paper is to present the different approaches used so far for the numerical approximation of nematic liquid crystal flows. Section 2 begins with a description of the function spaces which we will draw on throughout this work. We then describe the differential approaches commented above. In Section 3, we will analyze the advantages and disadvantages of existing finite element approximations in the literature for the Ginzburg-Landau penalization in terms of stability and convergence. We will make some remarks on the efficiency of the algorithms and linearization techniques in case of nonlinear schemes. Next, in Section 4, we review the numerical schemes designed for the saddle-point version. So far, the numerical analysis of these schemes is an open problem, mainly because the associated inf-sup condition for the director Lagrange multiplier is still not well-understood. 2 Problem statement 2.1 Some function spaces We will assume the following notation throughout this paper. Let Ω ⊂IRd, with d= 2 or 3, be a Lebesgue-measurable domain and let 1 ≤p≤ ∞. We denote by Lp(Ω) the space of all Lesbegue-measurable real-valued functions, f: Ω →IR, being pth-summable in Ω for p < ∞o or essentially bounded for p=∞, and by kfkLp(Ω) its norm. When p= 2, the L2(Ω) space is a 4 Hilbert space whose inner product is denoted by ·,·. Let α= (α1, α2, ..., αd)∈INdbe a multi-index with |α|=α1+α2+... +αd, and let ∂αbe the differential operator such that ∂α=∂ ∂x1α1...∂ ∂xdαd. For m≥0 and 1 ≤p≤ ∞, we define Wm,p(Ω) (see [1, 52]) to be the Sobolev space of m derivatives in Lp(Ω) whose norm is defined by kfkWm,p(Ω) = X |α|≤m k∂αfkp Lp(Ω)  1/p for 1 ≤p < ∞, kfkWm,p(Ω) = max |α|≤mk∂αfkL∞(Ω),for p=∞, where ∂αis understood in the distributional sense. In the particular case of p= 2, Wm,p(Ω) = Hm(Ω). Let C∞ c(Ω) be the space of infinitely differentiable functions with compact support in Ω. Then Wm,p 0(Ω) (analogously, Hm 0(Ω) when p= 2) is defined as the closure of C∞ c(Ω) in Wm,p(Ω). The dual spaces of Hs(Ω) and Hs 0(Ω) will be denoted by (Hs(Ω))′and H−s(Ω), respectively. For any space X, we shall denote the vector space Xdby its bold letter X. For example, (L2(Ω))dis denoted by L2(Ω), (Hm(Ω))dby Hm(Ω), etc. Consequently, in order to distinguish scalar-valued fields, such as the pressure p, from vector-valued fields, such as the velocity u, we denote them by roman letters and bold-face letters, respectively. For a real Banach space X,Lp(0, T ;X) denotes the space of X-valued functions fdefined on (0, T) such that k · kLp(0,T;X)=ZT 0 kfkp X1/p <∞.C1([0, T]; X) is the space of continuously differentiable X-valued functions in [0, T] such that supt∈[0,T]{kf(t)kX+kf′(t)kX}<∞. We will now introduce the function spaces in the context of the Navier-Stokes equations. Firstly, we define L2 0(Ω) = p:p∈L2(Ω),ZΩ p= 0, ϑ=nv∈C∞ c(Ω); ∇ · v= 0o. Then, let Hand Vbe the closure of ϑin L2(Ω) and H1 0(Ω), respectively, characterized by H={u∈L2(Ω); ∇ · u= 0,u·n= 0 on ∂Ω}, V={u∈H1(Ω); ∇ · u= 0,u=0on ∂Ω}. provided that Ω is Lipschitz [59]. 5 2.2 The Ericksen-Leslie equations Let Ω be a bounded open subset of Rd(d= 2 or 3) with boundary ∂Ω, and T > 0 the final time of observation. We will use the notation Q= Ω ×(0, T), Σ = ∂Ω×(0, T), and nthe unit outward normal to ∂Ω. The Ericksen-Leslie equations for the flow of a nematic liquid crystal can be written as                  ∂tu+ (u· ∇)u−ν∆u +∇p+λ∇ · ((∇d)t∇d) = 0in Q, ∇· u= 0 in Q, ∂td+ (u· ∇)d−γ∆d−γ|∇d|2d=0in Q, |d|= 1 in Q, (3) where u:Q→IRddenotes the solenoidal velocity field, p:Q→IR denotes the pressure, and d:Q→IRdrepresents the director field that describes the average molecular alignment. The parameter ν > 0 is a constant depending on the fluid viscosity, λ > 0 is an elastic constant and γ > 0 is a relaxation time constant. The expression (u·∇)uis the vector function whose ith component is Pd j=1 uj·∂jui,∇dis the gradient operator (∂jdi), (∇d)tdenotes its transpose and |d|=|d(x, t)|is the Euclidean norm in IRd. From the point of view of continuum mechanics, system (3) is essentially the simplest set of equations describing the motion of a nematic liquid crystal. System (3) was proposed by Lin in [41] from the macroscopic hydrodynamic theory of nematic liquid crystals established by Ericksen [22, 23] and Leslie [39, 40]. Since the original Ericksen-Leslie equations are mathematically untractable, further simplifications must be done in order to reduce the many reactive coupling terms between the fields dand uin the Oseen-Frank free energy functional. Although one could argue that system (3) is over-simplified, it keeps the core of the mathematical structure, such as strong nonlinearities and constraints, as well as the physical structure, such as the anisotropic effect of elasticity on the velocity vector field u. Thus, system (3) can be viewed as a good initial step towards the theoretical and numerical analysis of the original problem. Equation (3)1is the conservation of the linear momentum. The term λ∇ · ((∇d)t∇d) represents the anisotropic effect of the alignment on the fluid velocity. So, (3)1−2reduces to the classical Navier-Stokes equations for λ= 0, which indicates that the molecular centers of mass have no positional order; equation (3)2stands for the incompressibility of the fluid. Equation (3)3is the conservation of the angular momentum. Equation (3)4indicates that dis not a state variable, it only describes the orientation of the nematic liquid crystal molecules. As explained above, liquid crystals are in general sensitive to temperature. Herein, we only consider the thermally uncoupled model (3). To these equations we will add homogeneous 6 Dirichlet and Neumann boundary conditions for the velocity and director vector fields, respectively: u(x, t) = 0, ∂nd(x, t) = 0 on (x, t)∈Σ,(4) and the initial conditions u(x,0) = u0(x),d(x,0) = d0(x) on x∈Ω. (5) Here u0: Ω →IRdand d0: Ω →IRdare given functions such that ∇ · u0= 0 and |d0|= 1 in Ω. Pre-twist boundary conditions for the director dmodeled by non-homogeneous Dirichlet conditions as for the pixels could be straightforwardly considered but we have not included them here for the sake of clarity. A natural question that arises is: Does system (3) have an energy estimate? In other words, we want to know if system (3) is energetically dissipative. The next analysis is just done formally. Let us start observing that the elastic stress tensor can be re-written [42] as λ∇ · ((∇d)t∇d) = λ∇1 2|∇d|2−λ(∇d)t(−∆d),(6) where λ/2|∇d|2can be incorporated as part of the pressure. Then, multiply (3)1by uand (3)3 by λ(−∆d− |∇d|2d), and integrate over Ω. Provided that |d|= 1 can be proved previously in some way, it follows that λZΩ (∇d)t∆d·u−λZΩ (u· ∇)d·∆d= 0, λZΩ (u· ∇)d· |∇d|2d=λ 2ZΩ (|∇d|2u· ∇)|d|2= 0, and ZΩ ∂td· |∇d|2d=1 2ZΩ |∇d|2d dt|d|2= 0. As usual for the Navier-Stokes framework, the free-divergence constraint causes that pressure and convective term vanish together with the boundary condition for u. Therefore, one sees that problem (3) has the following energy law: d dt 1 2kuk2 L2(Ω) +λ 2k∇dk2 L2(Ω)+νk∇uk2 L2(Ω) +λγk∆d+|∇d|2dk2 L2(Ω) = 0.(7) Here we have used the homogeneous Dirichlet and Neumann boundary condition for the velocity and vector director field, respectively, to eliminate the boundary terms stemmed from integration by parts. Clearly, the sphere constraint over dhas been crucial to obtain the energy equality (7). This energy equality, which has been formally established, is called the first energy equality for (3). 7 It expresses the balance of energy in the system between the kinetic and elastic energies, i.e. the dependence between the linear and angular momentum equations along with the incompressibility and unit Euclidean norm constraint. Since any body force is considered, the rate of decay of the kinetic and elastic energy is dictated by the viscous term νk∇uk2 L2(Ω) and by the term λγk∆d+|∇d|2dk2 L2(Ω) which is the L2(Ω)-norm of the residual with respect to the steady equation (2). What one may now expect is that the constraint |d|= 1 would be a consequence of (3)1−3, such as the maximum principle for convection-diffusion equations, since it is an algebraic constraint. But this is not the case. First of all, we establish an equation for |d|. If we take the dot product of (3)3with d, this leads to 1 2 d dt|d(t)|2+1 2(u(t)· ∇)|d(t)|2−γ∆d(t)·d(t)−γ|∇d(t)|2|d(t)|2= 0. Next, observe that 1 2∆|d|2= ∆ d·d+|∇d|2.(8) Then, 1 2 d dt|d(t)|2+1 2(u(t)· ∇)|d(t)|2−γ 2∆|d(t)|2+γ|∇d(t)|2(1 − |d(t)|2) = 0. Equivalently, we have 1 2 d dt(|d(t)|2−1) + 1 2(u(t)· ∇)(|d(t)|2−1) −γ 2∆(|d(t)|2−1) −γ|∇d(t)|2(|d(t)|2−1) = 0.(9) Finally, if we multiply this equality by |d|2−1 and integrate over Ω, integration by parts yields that the convective term in this equality vanishes, since u∈V, and 1 4 d dt ZΩ (|d(t)|2−1)2+γ 2ZΩ (∇(|d(t)|2−1))2=γZΩ |∇d(t)|2(|d(t)|2−1)2.(10) In order to apply the Gronwall lemma, we would need to assume a regularity over dstronger than the one obtained from the energy estimate (7). To be more precise, ∇d∈L2(0, T;L∞(Ω)). Subsequently, we are not able to prove the maximum principle this way; and hence the sphere constraint must be imposed in system (3). The energy law (7) shows the proper functional spaces where a feasible definition of globalin-time weak solutions for (3) might be defined. Definition 1 A pair (u,d)is said to be a weak solution of (3)-(4)-(5) in (0, T)if: a) u∈L∞(0, T;H)∩L2(0, T;V), 8 d∈L∞(0, T;H1(Ω)), |d(x, t)|= 1,a.e. (x, t)∈Q. b) ∀φ∈C1([0, T]; V∩W1,∞(Ω)) such that φ(T) = 0, ZT 0 {−(u, ∂tφ) + ((u· ∇)u,φ) + ν(∇u,∇φ) −λ((∇d)t∇d,∇φ)dt = (u0,φ(0)). c) ∀ψ∈C1([0, T]; H1 0(Ω) ∩L∞(Ω)) such that ψ(T) = 0, ZT 0 {−(d, ∂tψ) + ((u· ∇)d,ψ) + γ(∇d,∇ψ) −γ(|∇d|2d,ψ)dt = (d0,ψ(0)). To our knowledge, the existence of a global-in-time weak solution for (3) still remains as an open problem. System (3) in its present formulation is overdetermined. There are four equations for three unknowns, and none of them seems to be dependent on the rest. Hence, the obtention of an energy estimates from (3) seems to be unaffordable. Of course, this is an important aspect when designing a numerical approximation for (3), because it could affect the robustness and uniqueness of the numerical method. As observed, equations (3)1(along with (3)2) and (3)3(along with (3)4) share a similar structure. But it is worth mentioning that (3)2is a linear differential constraint whereas (3)4 is a nonlinear algebraic constraint, which does not imply that the techniques in order to prove some results for one can be useful for the other. Of course, one expects that the constraints (3)2 and (3)4will be satisfied in some sense. With regard to (3)2, approximating solenoidal functions with conforming finite elements is difficult, in the sense that these finite-element basis functions are hard to construct and computationally inefficient. Thus, the treatment of the incompressibility (3)2, related to the pressure Lagrange multiplier, is enforced by means of velocity-pressure saddle-point methods which provide discrete divergence-free approximations and a well-known energy estimate. The corresponding velocity-pressure spaces must satisfy a discrete inf-sup condition, in order to get a well-posed discrete problem. This compatibility conditions between spaces can be circumvented via stabilized finite element methods (see e.g. [4]). Alternatively, the incompressibility condition can be penalized and the pressure unknown eliminated. This last approach does not require the computation of the pressure, but consistency is lost and the condition number blows up with the penalty, making the solution of the system too expensive for acceptable values of the penalty parameter (see e.g. [13]). 9 1. Stability. In numerical analysis, the terminology of stability is frequently associated with bounding the approximate solutions in certain norms, which are related to the energy of the model in question. Such estimates are called a priori or energy estimates. Moreover, the existence and uniqueness of approximate solutions is based on them. 2. Convergence. Two types of convergence results for approximate solutions can be considered. Compactness shows that the approximate solutions converge to a weak solutions under minimum regularity of the data, but no rate of convergence can be guaranteed. On the other hand, if a solution with more regularity than the one provided by the energy estimate is fixed, a priori error estimates can be performed, establishing rates of convergence. System (3)1−2becomes the classical Navier-Stokes equations plus an elastic stress tensor λ∇ · ((∇d)t∇d). Therefore, we cannot expect better results than those known for the Navier- Stokes equations (see [30, 59]). Since the numerical approximation of fluid flows is computationally expensive, efficient low-order approximations are favoured. Roughly speaking, the discrete velocity and pressure spaces are constructed by piecewise polynomials, globally continuous functions, satisfying the corresponding inf-sup condition. This sort of method provides algorithms which are easy to implement, well-conditioned with respect to the discretization parameters and allows to deal with complex geometries. An alternative to inf-sup stable approximations are the residual-based stabilization techniques (see e.g. [37]), that allow to choose the velocity and pressure spaces without the need to satisfy any compatibility condition and solves the singularly perturbed nature of the problem at hand for convection dominant flows. With regard to the approximation of (3)3−4, each method explained above will give rise to different variational formulations to be approximated by finite-dimensional spaces. In general, we will consider globally continuous piecewise polynomial functions, even though the velocity and director vector regularity are different. A notable exception is the method in [48], where C1 approximations are used to approximate the director vector; it was the first numerical scheme to deal with the approximation of (11). The Ginzburg-Landau method penalizes the constraint (7)4in the L2norm, say RΩ(|d|2− 1)2≤Cε2, from (12). Upon using this penalty method, an important choice is the size of the penalty parameter ε. A very high penalty number leads to ill-conditioned algebraic systems, when the off-diagonal blocks are multiplied by a large number. On the other hand, the rate of convergence is spoiled by ε. Sometimes, the penalization parameter depends on the mesh parameters. The Lagrange multiplier method introduces a new variable to be computed, which represents an increasement in the dimension of the resulting linear system. In any case, we will see that 16 the numerical approximations of the Ginzburg-Landau approach require auxiliary variables in order to prove an energy inequality, and in general the dimension of the resulting system is larger than the one of the saddle-point formulation. Furthermore, this method allows one to obtain numerical approximations of the original problem (without penalty) satisfying an energy inequality. On the other hand, the inf-sup condition must be satisfied in order to be well-posed (see [36]). As for the incompressibility condition, the constraint (11)4is satisfied in a discrete (weak) sense. Let us remark that, as far as we know, there are no finite element spaces capable to satisfy the restriction pointwise in general. E.g. we can straightforwardly prove that the only functions that belong to linear finite element spaces and satisfy the sphere constraint pointwise are constant functions. 2.3 Finite element approximation From now on, we assume that Ω is a bounded domain of IRd(d= 2 or 3) with a polygonal o polyhedral Lipschitz-continuous boundary and that there exists a family of triangulations {Th}h>0of Ω made up of triangles or quadrilaterals in two dimensions and tetrahedra or hexahedra in three dimensions, so that Ω = ∪K∈ThK. Further, any two elements K1, K2∈ Thsatisfy int(K1)∩int(K2) = ∅and K1∩K2is empty or a entire common vertex, face or side. For an arbitrary element K, we denote by hK>0 the diameter of K, with h= maxK∈ThhK, and by bKthe radius of the largest ball inscribed in K, with b= minK∈ThbK. The family of triangulations Thwill be assumed to be quasi-unform, i.e. there exists ρ > 0 such that b≥ρh for all K∈ Thand for all h > 0. The space of polynomials of degree less or equal to k > 0 in a finite element Kis denoted by Pk(K). The space of continuous piecewise polynomials is defined as Pk h=nvh∈C0(Ω) such that vh|K∈ Pk(K)∀K∈ Tho.(29) We also denote by Dk hthe space of piecewise polynomials of order no larger than kwithout C0 continuity. In particular, D0 his the space of piecewise constant functions. The notation of the finite element spaces to be used for approximating the primary variables are the following. For the velocity and pressure we let (Vh, Ph) to be two Lagrange finite element spaces associated with Th. Otherwise stated, the velocity finite element space Vhthat we consider is ((P1 h)d⊕ (Bh)d)∩H1 0(Ω), where Bh={vbsuch that vb|K∈ Pd+1(K), vb|∂K = 0,∀K∈ Th} is the space of bubbles in every element (see e.g. [10]). The pressure space Phis P1 h∩L2 0(Ω). 17 This velocity-pressure finite element pair is known to satisfy the discrete inf-sup condition kphkL2 0(Ω) ≤βsup v∈Vh\{0}qh,∇ · v kvkH1(Ω) ∀ph∈Ph,(30) for β > 0 uniform with respect to h. On the other hand, for the director field, we usually choose Dhto be again a Lagrange finite element space, i.e. (P1 h)d. Otherwise stated, these are the typical choices in the subsequent developments. In any case, the following results can be extended to any other finite element spaces, provided the required inf-sup conditions are satisfied. Throughout this work we will need to introduce some extra discrete spaces which will describe when necessary. We shall assume that Ω has the W2,r ×W1,r-elliptic regularity property for the Stokes problem. That is to say, for a prescribed f∈Lrwith r > 1, there exists a unique solution (u, p)∈W2,r(Ω) ×W1,r(Ω) of        −∆u+∇p=fin Ω, ∇ · u= 0 in Ω, u= 0 in ∂Ω, (31) which satisfies the following continuous dependence with respect to f: kukW2,r(Ω) +kpkW1,r(Ω) ≤CkfkLr(Ω).(32) On the other hand, we also assume that Ω holds a W2,r-elliptic regularity property for the Neumann problem; given g∈Lr(Ω) with RΩg= 0, there exists a unique solution d∈W2,r(Ω) with RΩd= 0 of (−∆d=gin Ω, ∂nd= 0 in ∂Ω,(33) which satisfies the following continuous dependence with respect to g: kdkW2,r(Ω) ≤CkgkLr(Ω).(34) For r= 2, properties (32) and (34) can be demonstrated for convex Ω [34, 38, 17, 18] or C1,1boundary [16, 55]. For r > 2, in two dimensions, properties (32) and (34) hold when Ω is convex, and rdepends on ∂Ω [34, 16]. However, in three dimensions, rdepends strongly on ∂Ω so that convexity is not suffice [16, 17]. An arbitrary triangulation of an arbitrary domain does not posses the C1,1regularity in general. Therefore, from the point of view of numerical analysis, we are limited to convex domains. Isoparametric elements are needed to generalize to domains with curved boundaries. 18 Usually, for the Navier-Stokes equation, the above regularity hypotheses are invoked when error estimates are stated, but we will see here that such a regularity is required even for convergence results by compactness. The initial data are u0∈Hand d0∈H1(Ω) with |d0|= 1 in Ω when the convergence of the algorithm is established by compactness or u0∈H2 0(Ω) ∩Vand d0∈W2,r(Ω) with |d0|= 1 in Ω when the convergence of the algorithm is established by error estimates. All the methods presented herein use a finite difference discretization in time. Let us therefore introduce some notation related to the time variable that we will use throughout the work. For simplicity, we suppose a uniform partition of [0, T ] into Npieces. So, the time step size is k=T/N and the time values (tn=nk)n=N n=0 . Let φn hbe a sequence in some Banach space X computed by some time-stepping scheme, i.e. the sequence φn hwill represent an approximation of φ(tn). We define δtφn+1 h=φn+1 h−φn h kand φn+1/2 h=φn+1 h+φn h 2. Let φr h,k and φl h,k be the piecewise constant interpolation taking the value φn+1 hor φn hon (tn, tn+1], respectively. Moreover, we define the piecewise linear interpolation φh,k ∈C0([0, T]; X) such that φh,k(tn) = φn h, that is, φh,k(t) = t−tj kφn+1 h+tj+1 −t kφn h∀t∈[tn, tn+1]. Let us introduce some short-hand notation, in order to simplify the writing of the different schemes: a(u,v) = (∇u,∇v) for all u,v∈H1(Ω), bu(u, q) = (∇ · u, q) for all u∈H1(Ω), q ∈L2(Ω), bd(q, d,¯ d) = hq, d·¯ difor all q∈(H1(Ω))′,d,¯ d∈H1(Ω) ∩L∞(Ω), c(w,u,v) = h(w· ∇)u·vifor all w,u,v∈H1(Ω). It is clear that c(w,u,u) = 0 for all w∈V. At the discrete level, the approximate velocity does not satisfy the incompressibility condition in a pointwise sense. Thus, in order to keep the skew-symmetry of the trilinear term, the stabilizing term 1 2h(∇ · u)w,viis added to the convective term. Then the trilinear form reads: ec(u,w,v) = h(u· ∇)w,vi+1 2h(∇ · u)w,vi,for all w,u,v∈H1(Ω), for which ec(w,u,u) = 0 holds for all u,w∈H1 0(Ω). 19 3 On the approximation of the Ginzburg-Landau problem (11) 3.1 Direct approximation for the Ginzburg-Landau problem 3.1.1 A H2-conforming approximation The first authors who dealt with the approximation of (11) for two-dimensional domains were Liu and Walkington in [48]. They focused their work in the obtention of error estimates, but they did not avoid the dependence on the penalty parameter ε. The best result for a Ginzburg- Landau-type equation is proved in [26], where the dependence of εis of polynomial order. Hence, an interesting open problem is to obtain error estimates independent of the penalty parameter ε, or at least a polynomial dependence. The starting point of the scheme lies in the following weak reformulation of (11), which is obtained by taking the L2inner product of ¯ u∈H1 0(Ω) with the linear momentum equation (11)3and the L2inner product of −∆¯ d(for ¯ d∈H2(Ω)) with the angular momentum equation (11)2. After using integration by parts, the problem reads as: find (u(t), p(t),d(t)) ∈H1 0(Ω) × L2 0(Ω) ×H2(Ω) satisfying                  (∂tu,v) + c(u,u,v) + νa(u,v) −bu(p, ¯ u) + λc(d,v,∆d) = 0, bu(¯p, u) = 0, (∂t∇d,∇¯ d)−c(u,d,∆¯ d) +γ(∆d−fε(d),∆¯ d) = 0, (35) for all (¯ u,¯p, ¯ d)∈H1 0(Ω) ×L2 0(Ω) ×H2(Ω). Observe that the elastic stress tensor has been written as in (21), where the potential term λ 2∇(|∇d|2) has been absorbed into the definition of the pressure, getting a modified pressure p∼p+λ 2|∇d|2. In [48], Fεwas truncated to have quadratic growth outside of the unit ball. In order to obtain a conforming approximation of (35) by a Galerkin finite element approximation, one should use a Dhspace based on H2conforming finite elements for the third equation in system (35). Thus, a general form for Dhis Dh={¯ dh∈C1(Ω) : ¯ d|T∈ Pd kfor all T∈ Th}. In two-dimensional domains, examples of C1finite element spaces are the bicubic Hermite elements or Argyris elements. However, for the first and second equation, one can use a classical pair (Vh, Ph) satisfying the inf-sup condition (30), as the one introduced above. In particular, Vh=P2 hand Ph=P1 hwere considered in [48]. With respect to the time integration, e.g. a fully implicit time stepping scheme is used in [48]. All these considerations lead to the following 20 scheme: given (un h, pn h,dn h), find (un+1 h, pn+1 h,dn+1 h)∈Vh×Ph×Dhsatisfying (δtun+1 h,¯ uh) + ec(un h,un+1 h,¯ uh) + ν a(un+1,¯ uh) −bu(pn+1 h,¯ uh) + λ c(¯ uh,dn+1 h,∆dn+1 h) = 0, bu(¯ph,un+1 h) = 0, (∇δtdn+1 h,∇¯ dh)−c(un+1 h,dn+1 h,∆¯ dh) +γ(∆dn+1 h−fε(dn+1 h),∆¯ dh) = 0, (36) for all (¯ uh,¯ph,¯ dh)∈Vh×Qh×Dh. As pointed out in [8], existence of a discrete solution (un+1 h, pn+1 h,dn+1 h) for (36) implies the relation between the time and penalty parameter k= O(e−1/ε2). This is due to the lack of an energy estimate independent of ε. In [48], the following discrete energy estimate for (36) was derived 1 2kun+1 hk2 L2(Ω) +λ 2k∇dn+1 hk2+k n X j=0 (νk∇uj+1 hk2+γλk∆dj+1 hk2) ≤C(1 2ku0 hk2+λ 2k∇d0 hk2) exp((n+ 1) k/ε2) for all n. (37) In order to obtain bounds for en u:= un h−u(tn) and en d:= dn h−d(tn), we must assume some sort of regularity for the solution of problem (11): u∈C(0, T;H2 0(Ω)), ∂tu∈L2(0, T;L2(Ω)), ∂ttu∈L2(0, T;H−1(Ω)) p∈C(0, T;H1(Ω)) and d∈C(0, T;W2,4(Ω)), ∂td∈L2(0, T ;H1(Ω)), ∂ttd∈L2(0, T;L2(Ω)). In addition, ∆d∈L2(0, T;H1(Ω)). In [48] it is shown the following error estimates: ken uk2 L2(Ω) +ken dk2 H1(Ω) +k n X m=0 kem uk2 H1(Ω) +kem dk2 H2(Ω)≤C(k2+h2) for all n, where the constant C > 0 depends on the penalty parameter εand, obviously, on the regularity of the solution. To obtain such error estimates, the differential form of (11) at t=tn+1 was used with the corresponding consistency errors. This approach needs the velocity vector u to satisfy an extra compatibility condition established in [35] on the data at t= 0 to get ∂ttu∈L2(0, T;H−1(Ω)). Moreover, the regularity for the director vector drequires the Neumann problem (33) holding the regularity (34) for r= 4, which leads to a restriction on the angles of the boundary ∂Ω. 21 Clearly, C1-finite elements provide high order approximations because they consider a huge degree of approximations. The Argyris element is constructed with polynomials of degree less than or equal to 5, and it has 21 degree of freedom per triangle for a scalar problem in dimension 2, whereas the bicubic Hermite element is obtained by products of polynomials of degree less than or equal to 3, involving 16 degree of freedom for the same case. The huge number of unknowns per element for a vectorial problem in two and three dimensions make this C1approximation extremely intensive in terms of computational cost. Furthermore, these approximations involve the derivatives of the unknowns, complicating their implementation. Although this scheme has demonstrated that can capture the behavior of singularities, it is not appropriate in terms of computation. The use of low order Lagrangian element is favoured for computational efficiency reasons. 3.1.2 A H1-conforming approximation The energy estimate (12) can be writen without the need to have d∈H2(Ω): d dt1 2kuk2 L2(Ω) +λ 2k∇dk2 L2(Ω) +λZΩ Fε(d) +νk∇uk2 L2(Ω) +λ γk∂td+ (u· ∇)dk2 L2(Ω) = 0, (38) since −∆d+fε(d), the equation for critical points of the energy RΩ 1 2|∇d|2+Fε(d), has been replaced by the material derivative ∂td+(u·∇)d. Lin and Liu presented in [44] one of the simplest time-stepping schemes for the two-dimensional Ginzburg-Landau problem (11), in which space is discretized by H1-conforming finite elements and time is discretized implicitly with respect to the linear terms and semi-implicitly with respect to the nonlinear terms, except for the anisotropic stress tensor that is fully explicit; so, the penalty term is discretized semi-implicitly. Thus, the resulting scheme reads (δtun+1 h,¯ uh) + ec(un h,un+1 h,¯ uh) + νa(un+1,¯ uh) −bu(pn+1 h,¯ uh)−λ((∇dn h)t∇dn h,∇¯ uh) = 0, bu(¯ph,un+1 h) = 0, (δtdn+1 h,¯ dh) + c(un h,dn+1 h,¯ dh) +γa(dn+1 h,¯ dh) + γ ε2((|dn h|2−1)dn+1 h,¯ dh) = 0, (39) for all (¯ uh,¯ph,¯ dh)∈Vh×Ph×Dh. In [44], the discrete spaces were Vh=P2 h,Ph=P1 hand Dh=P2 h. Since no energy estimates can be proved independent of the mesh and penalty parameter, the unique solvability of (39) is conditional, that is, there exist a polynomial relation k=R(ε, h) for which existence and uniqueness of discrete solutions for the scheme (39) may be established. 22 Obviously, this scheme reduces significantly the computational cost, allowing larger scale numerical simulations. The most important feature is that, at each time step, one only needs to solve a sequence of two decoupled linear problems for the velocity-pressure pair and the director field, separately. However, the authors were not able to derive the discrete analog to the energy estimate (38), which is basic for unconditional stability and also important with regard to the error analysis, if one wants to prove ken uk2 L2(Ω) +ken dk2 H1(Ω) +k n X m=0 kem uk2 H1(Ω) +kem dk2 H2(Ω)≤C(k2+h2) for all n. A second algorithm based on the finite element method of characteristics was presented in [44]. Fixed (un h, pn h,dn h), find (un+1 h, pn+1 h,dn+1 h)∈Vh×Ph×Dhsuch that (δtun+1 h(Xn h),¯ uh) + νa(un+1,¯ uh)−bu(pn+1 h,uh)−λ(∇dn h)t∇dn h,∇¯ uh= 0, bu(¯ph,un+1 h) = 0, (δtdn+1 h(Xn+1 h),¯ dh) + γa(dn+1 h,¯ dh) + γ ε2((|dn h|2−1)dn+1 h,¯ dh) = 0,(40) where we have denoted by δtun+1 h(Xn h) = un+1 h−un h(Xn h) kand δtdn+1 h(Xn+1 h) = dn+1 h−dn h(Xn+1 h) k, and Xn h:= x+kun h(x) and Xn+1 h:= x+kun+1 h(x), which is a backward Euler time discretization of the characteristic system X(x, tn+1;tn) = x−Ztn+1 tn u(X(x, tn+1;t), t)dt. (41) Note that the convective velocity for dn+1 his updated, due to the sequential feature of the method. The discrete spaces (Vh, Ph,Dh) are considered as above. The finite element method of characteristics for equations (11)1−2leads to linear algebraic problems with time-independent matrices, reducing the computational; we can decompose (at least approximately) the associated matrix just once at the beginning of computation. However, the same method does not avoid to compute for (11)3and (11)4, that lead to time-dependent matrices. For this reason, they proposed an iterative method for (40) in such a way the matrices do not change at each iteration. Given dn+1,j h, find dn+1,j h∈Dhsuch that dn+1,j+1 h−dn h(Xn+1 h) k,¯ dh!+a(dn+1,j+1 h,¯ dh) + γ ε2(|dn h|2−1)dn+1,j h,¯ dh= 0. It is proved the convergence of dn+1,j hto dn+1 hin the L2norm, but it requires the relation γk/ε2<1 between the time and penalty parameter. In any case, for large scale problems and 23 reasonably large time step sizes (e.g. ten times the explicit time step size), the assembling of the linear matrix is almost negligible compared to the solver time, and so, it is not clear the computational gain of this approach. On top of that, the iterations introduced over the d system are simple Richardson iterations, and so, it is expected to have a worse convergence than a Krylov-based solver. 3.2 Mixed Methods for the Ginzburg-Landau problem 3.2.1 Using w=∇das an auxiliary variable In a second work [49], Liu and Walkington avoided using Hermite finite elements for the approximation of the director equation in problem (11). The key idea is to introduce the auxiliary variable w=∇dwhich allows one to formulate the director equation in the framework of mixed methods. Then, (35) can be written as finding (u(t), p(t),d(t),w(t)) ∈H1 0(Ω)×L2 0(Ω)×H1(Ω)× H(div, Ω) such that                        (ut,v) + c(u,u,v) + νa(u,v) +bu(p, v) + λc(u,d,∇ · w) = 0, bu(p, u) = 0, (∂tw,¯ w)−c(u,d,∇ · ¯ w) +γ(∇ · w,∇ · ¯ w)−γ(fε(d),∇ · ¯ w) = 0, a(d,¯ d) + (∇ · w,¯ d) = 0, (42) for all (¯ v,¯p, ¯ d,¯ w)∈H1 0(Ω) ×L2 0(Ω) ×H1(Ω) ×H(div, Ω). Note that ∆d=∇ · w, which will be the extra equation to compute d, and the elastic stress tensor is written as (∇d)t∇ · wrather than wt∇·w; the convective term in the director equation remains the same (u·∇)d= (∇d)tw, rather than wtu. This is an important issue when obtaining error estimates, because one would need to establish high order norms for w, e.g. an L4norm, as used in [49]. However, this extra regularity is easier to obtain from the elliptic equation (42)4. Let us stress that whas d2 components, where dis the dimension space. A fully discrete scheme to approximate (42) is used in [49], which is implicit with respect to time: given (un h,dn h), compute (un+1 h, pn+1 h,wn+1 h,dn+1 h)∈(Vh, Ph,Wh,Dh) such that (δtun+1 h,¯ uh) + ec(un+1 h,un+1 h,¯ uh) + νa(un+1,¯ uh) −bu(pn+1 h,¯ uh) + λc(¯ uh,∇dn+1 h,∇ · wn+1 h) = 0, bu(¯ph,un+1 h) = 0, (δtwn+1 h,¯ wh)−c(un+1 h,dn+1 h,∇ · ¯ wh) +γ(∇ · wn+1 h−fε(dn+1 h),∇ · ¯ wh) = 0, a(dn+1 h,¯ dh)−(∇ · wn+1 h,¯ dh) = 0, (43) 24 for all (¯ uh,¯ph,¯ wh,¯ dh)∈Vh×Ph×Wh×Dh, where Whis an finite element approximation of H(Ω, div), e.g. the Raviart-Thomas [60] or the BDFM finite element [10]. The introduction of wallows to compute the director field and its first derivatives in an independent manner to keep an energy estimate similar to (37). However, the price to pay in terms of CPU cost is too large to be acceptable. The auxiliary unknown introduces four new components to be computed in dimension two, and nine in three dimensions, more than the number of unknowns of the original problem! Moreover, such a scheme is nonlinear, and some sort of iterative method must be performed at every time step. The solvability of scheme (43) by Brouwer fixed point argument is subject to satisfying k=O(e−1/ε2), which requires a very small time step for moderate values of ε. That restriction is a consequence of the energy estimate that scheme (43) provides, which is the same as for scheme (36), but expressed in terms of was follows: 1 2kun+1 hk2 L2(Ω) +λ 2kwn+1 hk2 L2(Ω) +k n X m=0 νk∇um+1 hk2 L2(Ω) +λ 2γk∇ · wm+1 hk2 L2(Ω) ≤Cexp((n+ 1) k/ε2)1 2ku0 hk2 L2(Ω) +λ 2kw0 hk2 L2(Ω)for all n, where the constant Cis a constant independent of ε. This bound blows up exponentially with ε−2. The hypotheses of regularity for a solution (u, p, d) of problem (11) are: u∈C(0, T;H2 0(Ω)), ∂tu∈L2(0, T;L2(Ω)), ∂ttu∈L2(0, T;H−1(Ω)), p∈C(0, T;H1(Ω)) and d∈C(0, T;W2,4(Ω)), ∂td∈L2(0, T ;H1(Ω)), ∂ttd∈L2(0, T;L2(Ω)). In addition, w∈C(0, T;H2(Ω)).Then the following error estimates hold: kun h−u(tn)k2 L2(Ω) +kwn h− ∇d(tn)k2 L2(Ω) +k N X n=0 kun h−u(tn)k2 H1(Ω) +k∇ · wn h−∆d(tn)kL2(Ω)≤C(k2+h2) for all n. As for scheme (36), the regularity required for ∂ttu∈L2(0, T;H−1(Ω)) needs a compatibility condition at t= 0 to be satisfied (see [35]). We would like to mention that this regularity assumption about the director vector to be approximated requires the domain Ω to hold (34) for the Neumann problem (33) with r= 4; therefore a restriction on the angles of the boundary must be imposed. 25 Proof: It suffices to prove (49) for all n= 0,...,N−1. For this, we argue by induction on n. Let us define Cd=K1+K2+K3with Kithe bounds for the initial data. Then, in particular, (u0 h,d0 h) satisfies the hypothesis of Lemma 3 for n= 0: ku0 hk2 L2(Ω) +λk∇d0 hk2 L2(Ω) + 2λRΩFε(d0 h)≤Cd, then (49) holds for n= 0 (and, in particular, ku1 hk2 L2(Ω) +λk∇d1 hk2 L2(Ω) + 2λRΩFε(d1 h)≤Cd). Now, we assume that (us h,ds h) holds (49) for s= 1, ..., n −1. Adding (49) for s= 1, ..., n −1, one has kun hk2 L2(Ω) +λk∇dn hk2 L2(Ω) + 2λZΩ Fε(dn h)≤ ku0 hk2 L2(Ω) +λk∇d0 hk2 L2(Ω) + 2λRΩFε(d0 h) ≤K1+K2+K3=Cd, which implies from Lemma 3 that (49) holds for n. The convergence for scheme (44) is demonstrated by two different ways: compactness and error estimates. The compactness for the discrete velocity is attained by estimating a discrete fractional time estimate for the velocity and director field. The line of argument is as follows. Add (44)1from n=mto m+r−1 and take ¯ u=um+r h−um has a test function. Then, use adequately equation (44)2and add from m= 1 to N−rto get k N−r X m=1 kum+r h−um hk2 L2(Ω) =−ν k2 N−r X m=1 m+r X n=m a(un+1 h,um+r h−um h) −k2 N−r X m=1 m+r X n=m ch(un h,un+1 h,um+r h−um h) −λ k2 N−r X m=1 m+r X n=m c(um+r h−um h,dn h,wn+1 h). Let us now focus on how to control the last term. A discrete integration by parts and Sobolev’s inequality leads to λk2 N−r X m=1 m+r X n=m c(um+r h−um h,dn h,wn+1 h) ≤C(k r)1/2 N−1 X n=1 kkwn+1 hk2 L2(Ω)!1/2 N−1 X n=1 kk∇ dn hk2 L3(Ω)!1/2 N−r X m=1 kk∇(um+r h−um h)k2 L2(Ω)!1/2 ≤C(k r)1/2.(61) Obviously, we need the control kPN−1 n=1 kdn hk2 W1,3(Ω) ≤Cto obtain (61), which takes advantage of kPN−1 n=1 kwn hk2 L2(Ω) ≤C, proved in [27] by assuming ∂Ω to be only Lipschitz. Therefore, we arrive at N−r X m=0 kum h−um+r hk2 L2(Ω) dt ≤C(rk)1/2for all r= 1, ..., N −1. 32 A similar bound for the discrete director can be obtained: N−r X m=0 kdm+r h−dm hk2 H1(Ω) dt ≤C(rk)1/2for all r= 1, ..., N −1. The expressions above, in term of time interpolations associated with the sequences {un h}and {dn h}, are written as ZT−δ 0 kur h,k(t+δ)−ur h,k(t)k2 L2(Ω) dt ≤C δ1/2for all δ∈(0, T ),(62) ZT−δ 0 kdr hk(t+δ)−dr hk(t)k2 H1(Ω) dt ≤C δ1/2for all δ∈(0, T).(63) Next, a compactness result from [58] provides that the sequences {dr h,k}and {ur h,k}are compact in L2(0, T;L2(Ω)), which is extensible to {dl h,k}and {ul h,k}. This strong convergence for {dr h,k}is not enough to pass to the limit in the elastic stress tensor of equation (44)1. Testing equation (44)4with ¯ d=dn+1 hone can prove the strong convergence of {∇dr h,k}and {∇dl h,k}towards ∇din L2(0, T;L2(Ω)) as obtained in [27]. Therefore, one gets the existence of a global-in-time weak solution to (11) (see Definition 2). On the other hand, error estimates are obtained in [27]. Indeed, if the solution (u, p, d) to (11) has the following regularity: (u,d,w)∈L2(H2(Ω) ×W2,3(Ω) ×H1(Ω)) ,u∈L∞(W1,3(Ω) ∩L∞(Ω)) , p ∈L2(H1(Ω)), (∂tu, ∂td, ∂tw)∈L2(H1(Ω) ×H2(Ω) ×W1,6/5(Ω)), then scheme (44) satisfies the error estimates: ken uk2 L2(Ω) +k∇en dk2 L2(Ω) +k n+1 X m=1 k∇em uk2 L2(Ω) +kem wk2 L2(Ω)≤C(h2+k2) for all n. where we recall that en u:= un h−u(tn), en d:= un h−d(tn), and en w:= wn h−w(tn). These error estimates are derived by using an integral formulation for problem (11) which avoids to assume nonlocal compatibility conditions on the initial data [35]. Instead, the integral formulation requires more regularity for the time derivatives of the solution to be approximated. Three iterative methods to decouple (un+1 h, pn+1 h) from (dn+1 h,wn+1 h) at each time step are also obtained for scheme (44) in [27]. For any n≥0, given wn h,dn h,un hand pn h, these schemes compute wn+1 h,dn+1 h,un+1 hand pn+1 has follows: 1. (a) Given w0=wn h, if n > 1, or w0=−PWh(∆d0), if n= 1, where PWhis the L2 orthogonal projection onto Wh. 33 (b) For i≥1, known wi−1, compute (ui, pi)∈Vh×Phsuch that 1 k(ui−un h,¯ u) + νa(∇ui,∇v)−bu(pi,∇ · v) +ch(un h,ui,¯ uh) = c(¯ uh,dn hwi−1), bu(¯ph,ui) = 0. (c) Next, compute (wi,di)∈Wh×Dhthe solution of: 1 k(di−dn h,¯ wh) + γ(wi,wh) = −c(ui,dn h,¯ wh)−γ(e fε(dn h),¯ wh), a(di,¯ dh)−(wi,¯ wh) = 0. 2. (a) Let u0=un h. (b) Known ui−1, compute (wi,di)∈Wh×Dh, solution of: 1 k(di−dn h,¯ wh) + γ(wi,¯ wh) = −c(ui−1,dn h,¯ wh)−γ(e fε(dn h),¯ wh), (di,¯ dh)−(wi,¯ dh) = 0. (c) Next, compute (ui, pi)∈Vh×Phsuch that 1 k(ui−un h,¯ uh) + νa(ui,¯ uh)−bu(pi,¯ uh) +ch(un h,ui,¯ uh) = c(¯ uh,dn h)t,wi), bu(¯ph,ui) = 0. 3. (a) Let u0=un hand w0=wn h, if n > 0, or w0=−PWh(∆d0), if n= 1. (b) Known ui−1and wi−1, compute in a parallel way •(wi,di)∈Wh×Dhsuch that 1 k(di−dn h,¯ wh) + γ(wi,¯ dh) = −c(ui−1,dn h,¯ wh)−γ(e f(dn h),¯ wh), a(di,¯ dh)−(wi,dh) = 0, •and (ui, pi)∈Vh×Phsuch that 1 k(ui−un h,¯ uh) + νa(ui,¯ uh)−bu(pi,uh) =−ch(un h,ui−1,¯ uh) + c(¯ uh,dn h,wi−1), bu(¯phui) = 0. The convergence of these iterative methods is established under the condition h2≤Ck, where C > 0 is a constant depending on ε. 34 3.2.3 Using w=−∆d+fε(d)as an auxiliary variable So far, the finite element schemes that have been presented to approximate a solution to the Ericksen-Leslie equations (3) by means of the Ginzburg-Landau equations (11) do not preserve a discrete version of the energy law (12). In fact, the only exception is the novel modification of (44), i.e. system (47). This is basically due to the fact that the nonlinear function fε(d) does not belong to the discrete space Dh. In order to get such an energy law one must write the elastic tensor λ∇ · ((∇d)t∇d) in terms of the critical point equation −∆d+fε(d) like λ∇ · ((∇d)t∇d) = λ∇1 2|∇d|2+Fε(d)−λ(∇d)t(−∆d+fε(d)),(64) as in [8, 33]. Then, if one defines the variable w=−∆d+fε(d), which represents the Euler- Lagrange equation related to the Ginzburg-Landau free energy functional RΩ 1 2|∇d|2+Fε(d), model (11) can be reformulated as              ∂tu+u· ∇u−ν∆u+∇p−λ(∇d)tw=0in Q, ∇· u= 0 in Q, |d| ≤ 1, ∂td+u· ∇d+γw=0in Q, −∆d+fε(d)−w=0in Q, (65) where the pressure pis modified by the potential function p+λ 2|∇d|2+λFε(d) (which is called again pfor simplicity). The vector spaces where the weak formulation of problem (65) is well-posed are as follows: find (u(t), p(t),d(t),w(t)) ∈H1 0(Ω) ×L2 0(Ω) ×H1(Ω) ×L2(Ω) such that (∂tu,¯ u) + c(u,u,¯ u) + ν a(u,¯ u) −bu(p, ¯ u)−λ c(¯ u,d,w) = 0, bu(¯p, u) = 0, (∂td,¯ w) + c(u,d,¯ w) + γ(w,¯ w) = 0, a(d,¯ d) + (fε(d),¯ d)−(w,¯ d) = 0, (66) for all (¯ u,¯p, ¯ w,¯ d)∈H1 0(Ω) ×L2 0(Ω) ×L2(Ω) ×H1(Ω). Two finite-element Euler time-stepping schemes [8, 33] have been developed to approximate (66), being both implicit for the linear terms and semi-implicit for the nonlinear ones. The main difference among them lies in the way of treating the time integration of fε. Becker, Feng and Prohl proposed in [8] a fully implicit approximation while Guill´en-Gonz´alez and Guti´errez- Santacreu suggested a fully explicit one. Thus, the scheme developed in [8] is expressed as follows. Let (un h,dn h)∈(Vh,Dh) be given, 35 then find the solution (un+1 h, pn+1 h,dn+1 h,wn+1 h)∈Vh×Ph×Dh×Whof the nonlinear system:                  (δtun+1 h,¯ uh) + ec(un h,un+1 h,¯ uh) + νa(un+1 h,¯ uh) −bu(pn+1 h,¯ uh)−λc(¯ uh,dn h,wn+1 h) = 0, bu(¯ph,un+1 h) = 0, (δtdn+1 h,¯ wh) + c(un+1 h,dn h,¯ wh) + γ(wn+1 h,¯ wh) = 0, a(dn+1 h,¯ dh) + (fε(dn+1 h,dn h),¯ dh)h−(wn+1 h,¯ dh) = 0, (67) for all (¯ uh,¯ph,¯ wh,¯ dh)∈Vh×Ph×Wh×Dh, where fε(dn+1 h,dn h) = 1 2ε2|dn+1 h|2dn+1 h−dn h. Here, a discrete inner product is used (·,·)h, which is defined as follows. Let {φa:a∈Nh}denote the nodal basis associated the set of all nodes Nh={al}l∈Lof Th. Thus, the nodal interpolation operator IDh:C(Ω) →Dhis such that IDhψ:= Pa∈Nhψ(a)φa. Then the discrete inner product (·,·)his defined in the following way: for all Ψ,ψ∈C(Ω), one has (Ψ,ψ)h:= ZΩ Ih(Ψ·ψ) = X a∈Nh Ψ(a)·ψ(a)ZΩ φa. This discrete inner-product applied to the potential term (fε(dn+1 h,dn h),¯ dh)hproduces a 1/ε2- diagonal lumped mass matrix. A Newton method is considered in [8] to linearize the problem. Scheme (67) provides a discrete energy law which mimics the continuous energy law (12): 1 2kun+1 hk2 L2(Ω) +λ 2k∇dn+1 hk2 L2(Ω) +λZΩ Fε(dn+1 h) + k n+1 X m=1 νk∇um hk2 L2(Ω) +λγkwm hk2 L2(Ω) ≤1 2ku0 hk2 L2(Ω) +λ 2k∇d0 hk2 L2(Ω) +ZΩ Fε(d0 h)for all n. (68) This scheme is unconditionally stable and convergent towards a measured-valued solution of (3). This convergence is attained in two steps; firstly, when the time and space discretization parameters go to zero, the convergence towards a weak solution of the penalized problem (65) is proved, and afterwards, when the penalty parameter εgoes to zero, one arrives at a measurevalued solution of problem (3), where the elastic tensor (∇d)t∇donly tends to a certain measure (see [43]). The scheme presented in [33] was designed for nonhomogeneous Dirichlet boundary conditions. This scheme should be redesigned appropriately in order to guarantee stability for homogeneous Neumann boundary conditions. The renewed scheme preserves the time discretization of the penalty function in a fully explicit way. Then, the scheme becomes: given (un h,dn h)∈(Vh,Dh), find (un+1 h, pn+1 h,dn+1 h,wn+1 h)∈Vh×Ph×Dh×Whsolving the finite element linear system:                  (δtun+1 h,¯ uh) + ec(un h,un+1 h,¯ uh) + νa(un+1 h,¯ uh) −bu(pn+1 h,¯ uh)−λc(¯ uh,dn h,wn+1 h) = 0, bu(¯ph,un+1 h) = 0, (δtdn+1 h,¯ wh) + c(un+1 h,dn h,¯ wh) + γ(wn+1 h,¯ wh) = 0, a(dn+1 h,¯ dh) + (fε(dn h),¯ dh)h−(wn+1 h,¯ dh) = 0, (69) for all (¯ uh,¯ph,¯ wh,¯ dh)∈Vh×Ph×Wh×Dh. Scheme (69) is conditionally stable in the terms of scheme (47) (assuming (S)) and convergent (assuming (C)) to a measure-valued solution, but this time the convergence is attained by making the mesh parameter (h, k) and the penalty parameter εgo to zero at the same time, what could be done for scheme (67) as well. Since the approximation of fεis explicit, the linear algebraic problem (69) does not depend on ε, which assures that the conditional number of the system at every time step is not affected by ε. However, the stability condition (S) implies that the time step kmust be quite small if the size of εis proportional to the space parameter h, but numerical experiences in [33] have demonstrated to be optimal. The convergence of schemes (67) and (69) is obtained by means of compactness results. For the discrete director vector, it is easy to check that k N−1 X n=0  δtdn+1 h 2 L3/2(Ω) ≤C, by choosing as a test function ¯ wh=PWh¯ win (67) or (69); we have also used the fact that w∈ L3(Ω) and the L3(Ω)-stability of the L2projection operator PWh(see [19]). As a consequence of the energy estimates and a compactness result in [58], one gets the compactness of the sequence {dh,k,ε}in Lq(0, T;Lr(Ω)) with 1 ≤r < 6 and 1 ≤q < ∞, where dh,k,ε is the linear piecewise continuous function such that dh,k,ε(tn) = dn h. We think the way how the compactness for the discrete velocity in L2(0, T;L2(Ω)) is proved in [8] is not clear. From (68), the sequence of discrete velocities is bounded in L∞(0, T;L2(Ω))∩ L2(0, T;H1 0(Ω)). Afterwards, by a duality argument, the discrete time derivative for the velocity is bounded in the dual space of V∩H2(Ω). Then the Aubin-Lions compactness lemma is used. To apply this compactness result, one needs the embeddings H1 0(Ω) compact 7→ L2(Ω) ֒→ (V∩H2(Ω))′, but the embedding from L2(Ω) into (V∩H2(Ω))′is not injective. Obtaining a compactness result for the discrete velocity turns out to be harder than for the Ginzburg- Landau problem (11) for a fixed ε, even for problems with a similar structure such as the density-dependent Navier-Stokes equations ([50], [31] and [32]). 37 Now we cannot prove the time fractional estimate N−r X m=0 kum h−um+r hk2 L2(Ω) dt ≤C(rk)1/2for all r= 1, ..., N −1 obtained in [27], since no control over kPN−1 n=1 kdn hk2 W1,3(Ω) ≤Cis available with Cbeing independent of ε. Let us sketch the way of getting compactness for the discrete velocity in an appropriate way. Let Xh={¯ uh∈Vh:bu(¯ph,¯ uh) = 0 ∀¯ph∈Ph} be the discrete divergence-free space associated with Vhand consider A−1 h:Vh→Xhthe inverse discrete Stokes operator defined as (∇A−1 huh,∇vh) = (uh,vh)∀vh∈Xh.(70) Notice that (70) is well-defined owing to the inf-sup condition (30). By multiplying (67)1by k2, summing for n=m, ..., m−1+r, setting ¯ uh=A−1 h(um+r h−um h) as a test function in (67)1, and summing for m= 0, ..., N −r, we get k N−r X m=0 |∇A−1 h(um+r h−um h)|2=−k2 N−r X m=0 m−1+r X n=m cun h,un+1 h, A−1 h(um+r h−um h) +ν k2 N−r X m=0 m−1+r X n=m∇un+1,∇A−1 h(um+r h−um h) +λ k2 N−r X m=0 m−1+r X n=m(∇dn h)twn+1 h, A−1 h(um+r h−um h) := J1+J2+J3. (71) Let us only focus on how to estimate J3 J3≤C k2 N−r X m=0 m−1+r X n=m |∇dn h||wn+1 h|kA−1 h(um+r h−um h)kL∞(Ω) ≤C k2 N−r X m=0 m−1+r X n=m |wn+1 h|kA−1 h(um+r h−um h)kL∞(Ω). Sobolev’s inequality shows that kA−1 h(um+r h−um h)kL∞(Ω) ≤CkA−1(um+r h−um h)kW1,r(Ω) with r > d,dbeing the dimension of Ω. The following bound kA−1 h(um+r h−um h)kW1,r(Ω) ≤CkA−1(um+r h−um h)kW1,r(Ω), 38 is proved in [28], under the regularity (32) with r > d for A−1the Stokes operator (31). Thus, applying Sobolev’s inequality, H2(Ω) ֒→W1,r(Ω), with r≤6, gives kA−1(um+r h−um h)kL∞(Ω) ≤CkA−1(um+r h−um h)kH2(Ω) ≤Ckum+r h−um hkL2(Ω),(72) where we have used the regularity result (32) for r= 2. Therefore, J3≤C k2 N−r X m=0 m−1+r X n=m kwn+1 hkL2(Ω)kum+r h−um hkL2(Ω) ≤k N−r X m=0 kum+r h−um hkL2(Ω) k m−1+r X n=m kwn+1 hk2 L2(Ω)!1/2 (r k)1/2≤C(r k)1/2. Finally, we conclude that k N−r X m=0 k∇A−1 h(um+r h−um h)k2 L2(Ω) ≤C(r k)1/2, which is equivalent to ZT−δ 0 kuh,k,ε(t+δ)−uh,k,ε(t)k2 X′ hdt ≤C δ1/2∀δ: 0 < δ < T. due to the fact that k∇A−1 huhkL2(Ω) and kuhkX′ hare equivalent norms. Note that the time fractional time estimate for the discrete velocity is bounded in the V′ h norm, which depends on the space parameter h; therefore we cannot apply the compactness results given by Simon in [58]. Then the idea will be to encounter a fractional time norm being independent of the mesh parameters. Consider the space V={u∈H1 0(Ω) : ∇ · u= 0}and the orthogonal projection Rh:Vh→Vsuch that ∇(Rhuh−uh),∇v= 0,∀v∈V. One knows from [33] that kRhuhkV′≤Ch|∇ · uh|+kuhkV′ hand ZT−δ 0 kRhuh,k,ε(t+δ)−Rhuh,k,ε(t)k2 V′dt ≤C δ1/2+C h. Finally, the compactness of {Rhul h,k,ε}in L2(0, T;L2(Ω)) follows by a perturbed compactness result due to Az´erad and Guill´en-Gonz´alez in [3]. To conclude with the strong convergence, one uses the external approximation from Xhto V, in order to prove that the sequence {ul h,k,ε} is compact in L2(0, T;L2(Ω)). To complete with the convergence we must pass to the limit. In [33], the authors needed to impose the additional hypothesis (C) in the process. Although, scheme (67) is unconditionally stable, it may need such a condition. An interesting issue is to prove compactness under the minimum assumptions over the boundary of Ω. At the continuous level, one would need the Neumann problem to have the regularity (34) for r= 2 when a semi-Galerkin method is used, while the discrete compactness result needs the Stokes problem to have the regularity (32) for r > d, with dthe space dimension. For the Ginzburg-Landau problem (11), it was attained in [27] for fixed ε. 39 3.2.4 A second-order scheme for the Ginzburg-Landau problem We next discuss the work of Lin, Liu and Zhang in [45]. They presented the first numerical work preserving an exact discrete energy law for (38), i.e. without introducing numerical diffusion provided by the time-stepping schemes. It relies on a modified Crank-Nicolson or midpoint scheme: given (un h,dn h)∈(Vh,Dh) , find (un+1 h,wn+1 h)∈Vh×Whand (pn+1 h,dn+1 h)∈Ph×Dh solving the finite element nonlinear system:                        (δtun+1 h,¯ uh) + ec(un+1/2 h,un+1/2 h,¯ uh) +νa(un+1/2 h,¯ uh)−bu(pn+1/2 h,¯ uh) −λ γc(¯ uh,dn+1/2 h, δtdn+1 h+ (un+1/2 h· ∇)dn+1/2 h) = 0, bu(¯ph,un+1/2 h) = 0, (δtdn+1 h,¯ dh) + c(un+1/2 h,dn+1/2 h,¯ dh) +γa(dn+1/2 h,¯ dh) + γ(fε(dn+1 h,dn h),¯ dh) = 0, (73) where fε(dn h,dn+1 h) = 1 ε2 (|dn+1 h|2−1) + (|dn h|2−1) 2dn+1/2 h. Clearly, this approximation of fε(d) has second order of accuracy in time. The energy law given by (73) takes the form 1 2kun+1 hk2 L2(Ω) +λ 2k∇dn+1 hk2 L2(Ω) +λZΩ Fε(dn+1 h) +ν k n+1 X m=1 k∇um+1/2 hk2 L2(Ω) +λ γk n+1 X m=1 kδtdm+1 h+ (um+1 2 h· ∇)dm+1 2 hk2 L2(Ω) =1 2ku0 hk2 L2(Ω) +λ 2k∇d0 hk2 L2(Ω) +ZΩ Fε(d0 h) for all n. The fact that no auxiliary variable is introduced is a clear benefit of this approach. In contract, since scheme (73) is nonlinear, one needs to perform nonlinear iterations. The authors proposed a Picard type linearization together with a block Gauss-Seidel solver and a fully explicit treatment of the nonlinear terms, in order to have a time-independent time matrix and use an exact LU or Cholesky decomposition only once. Let us define φs−1/2=φn h+φs−1 h 2and φs+1/2=φs h+φn h 2. Then, given (un h,dn h)∈(Vh,Dh), the methods reads as: 1. Known us−1, find (us, ps)∈Vh×Phsuch that (us h−un h k,¯ uh) + ec(us−1/2,us−1/2,¯ uh) +νa(us+1/2,¯ uh)−bu(ps+1/2,¯ uh) −λc(¯ uh,ds−1/2,ds−1 h−dn h k+ (us−1/2· ∇)ds−1/2) = 0, bu(¯ph,us+1/2) = 0. 40 2. Next, find ds∈Dhsuch that (ds h−dn h k,¯ dh) + c(us−1/2,ds−1/2,¯ dh) +a(ds+1/2,¯ dh) + (fε(dn h,ds−1 h),¯ dh) = 0. The use of direct solvers is restricted to very small problems. For real applications, iterative methods are the only choice, due to CPU and memory limitations. In those cases, this way to deal with coupling and nonlinearities is very rudimentary. More modern approaches to linearization and linear solvers, e.g. Newton-Krylov-type solvers, surely provide better results. There is no numerical analysis for scheme (73). The convergence in the sense of error estimates or compactness is a interesting question that is open. 4 On the approximation of the Ericksen-Leslie problem (3) This section is devoted to finite element schemes for the Ericksen-Leslie equations (3) that provide a discrete energy law. To obtain the energy law (7) we need the sphere constraint |d|= 1 to be fulfilled almost everywhere in Qwhich is difficult to achieve at the discrete level due to the own nature of Lagrange finite elements. Therefore, constructing stable numerical approximations directly for (3) results a more difficult task than for the Ginzburg-Landau problem (11). 4.1 A direct approximation The next scheme we will present is due to Becker, Feng, and Prohl in [8]. The numerical approximations are based on the ideas given in [7] for the unsteady harmonic map equation which utilizes the Galerkin method with Lagrange finite elements of order 1. The vector identity a×(b×c) = (a·c)b−(a·b)cfor all a,b,c∈IR3and |d|= 1 lead to the identity d×(d×∆d) = −|∇d|2d−∆d. Therefore, problem (3) becomes                  ∂tu+u· ∇u−ν∆u +∇p+λ∇ · ((∇d)t∇d) = 0in Q, ∇· u= 0 in Q, ∂td+u· ∇d+γd×(d×∆d) = 0in Q, |d|= 1 in Q. (74) From equation (74)3, it is not hard to prove that |d|= 1 holds almost everywhere in Q. Therefore, the set of equations (3)3−4is equivalent to equation (74)3. Based on that, we consider 41 [2, 6] is also performed. The idea is to consider a Picard linearization of (81) with (82) together with a normalization of the director field approximation. Let (un+1,k h, pn+1,k h,dn+1,k h, qn+1,k h,)∈Vh×Ph×Dh×Qhbe known. 1. Compute (˜ dn+1,k+1 h, qn+1,k+1 h)∈Dh×Qh, solution of (δt˜ dn+1,k+1 h,¯ dh) + c(un+1,k h,˜ dn+1,k+1 h,¯ dh) +γa(˜ dn+1,k+1 h,¯ dh) + γbd(qn+1,k+1 h,dn+1,k h,¯ dh) = 0, bd(¯qh,dn+1,k h,˜ dn+1,k+1 h) = h1,¯qhi. 2. Compute dn+1,k+1 h∈Dhas dn+1,k+1 h(a) = ˜ dn+1,k+1 h(a) |(˜ dn+1,k+1 h(a)| for all a∈ Nh. 3. Known dn+1,k+1 h∈Dh, compute (un+1,k+1 h, pn+1,k+1 h)∈Vh×Phby a Picard linearization: (δtun+1,k+1 h,¯ uh) + ec(un+1,k h,un+1,k+1 h,¯ uh) + νa(un+1,k+1 h,¯ uh) +bu(pn+1,k+1 h,vh) + λ γc(¯ uh,dn+1,k+1 h, ∂tdn+1,k+1 h+ (un+1,k+1 h· ∇)dn+1,k+1 h) = 0, bu(¯ph,un+1,k+1 h) = 0. Obviously, dn+1,k+1 hsatisfies the nonlinear constraint (87) since Step 2 enforces the sphere constraint on the nodes. With regard to CPU cost, the saddle-point structure needs eight degrees of freedom per node, only beated by schemes (39), (73), and (77) which do not introduced any auxiliary variable. The numerical analysis of schemes (81), (83), and (84) is still open. 5 Conclusions In this work, we have analyzed the existing numerical schemes in the literature to approximate the Ericksen-Leslie equations (3) by means of low-order finite elements. We have distinguished between implicit and semi-implicit methods. Existence of discrete solutions have been detailed for the implicit schemes, as well as their linearization. The convergence of these algorithms is presented by two different ways: compactness and error estimates. Clearly, schemes (36) and (43) have been designed to satisfy the energy law (17). This is the reason why we think that there is no way to adapt them to hold an energy estimate independent of the penalty parameter ε. The large number of degrees of freedom make them somehow impracticable for large scale simulations. 48 Probably, scheme (39) is the most efficient algorithm among the linear methods presented herein since it does not compute any extra auxiliary variable. In contrast, there is no mean of finding an energy estimate independent of ε. Schemes (47), (67) and (69), designed for approximating the Ginzburg-Landau equation, need to compute extra variables in order to keep a discrete version of the energy law (7). They introduce between dand d2additional degrees of freedom per node (dbeing the space dimension), which implies a too high computationally cost. In particular, although scheme (69) has such an energy estimate, there is no control over the auxiliary unknown w=−∆dindependent of ε, which could clearly deteriorate the resulting linear system. Another important drawback of schemes (47) and (69) is the fact that they are conditionally stable. The relation (45) is quite restrictive, in the sense that one needs a time step small for moderate values of εif ε=O(h). It is interesting to note that scheme (67) is unconditionally stable, but it requires the extra relation (46) to pass to the limit towards measure-valued solutions, as was pointed out in [33]. One of the main advantages of using the saddle-point approach is that it allows to approximate numerically both the Ericksen-Leslie problem and the Ginzburg-Landau problem by means of the same numerical approximation. Moreover, it allows to take εarbitrary small (or even zero) in comparison with the previous algorithms, since the condition number is independent of ε. However, the proof of the inf-sup condition that would allow to prove the well-posedness of the algorithm is still open. Numerical tests have shown an (H1)′stability in [5]. We would also mention that the nonlinearity in the equation for the director vector is quadratic while for other alternatives it is cubic. The best choice to approximate the Eriscksen-Leslie equations in terms of the number of degrees of freedom is scheme (77), since it does not introduce any extra variable. Nevertheless, the convergence is not clear even for two-dimensional problems due to the fact the cubic nonlinearity involves second derivatives. References [1] R. A. Adams. Sobolev Spaces. Academic Press, New York, 1975. [2] F. Alouges. A new algorithm for computing liquid crystal stable configurations: the harmonic mapping case. SIAM Journal on Numerical Analysis, 34 (5), 1708-1726, 1997. [3] P. Az´ erad, F. Guill´ en-Gonz´ alez.Mathematical justification of the hidrostatic approximation in the primitive equations of geophysical fluid dynamics. SIAM J. Math. Anal. 33 (4), 847-859, 2001. 49 [4] S. Badia and R. Codina.On a multiscale approach to the transient Stokes problem: Dynamic subscales and anisotropic space-time discretization. Applied Mathematics and Computation 207, 415-433, 2009. [5] S. Badia, F. Guill´ en-Gonz´ alez, J. V. Guti´ errez-Santacreu.Finite element approximation of nematic liquid crystal flows using a saddle-point structure. Submitted. [6] S. Bartels.Stability and convergence of finite-element approximation schemes for harmonic maps. SIAM Journal on Numerical Analysis, 43(1), 220-238, 2005. [7] S. Bartels, A. Prohl.Constraint preserving implicit finite element discretization of harmonic map heat flow into spheres. Math. Comp., 76, 1847-1859, 2007. [8] R. Becker, X. Feng, A. Prohl. Finite element approximations of the Ericksen-Leslie model for nematic liquid crystal flow. SIAM J. Numer. Anal. 46(4), 1704-1731, 2008. [9] S. Brenner, L. R. Scott.The Mathematical Theory of Finite Element Methods, TAM 15, Springer-Verlag, Berlin, 1994. [10] F. Brezzi and M. Fortin. Mixed and hybrid finite element methods, no. 15 in Computational Mathematics. Springer-Verlag (1991). [11] Y.M. Chen. The weak solutions to the evolution problems of harmonic maps. Math. Z. 201(1), 69-74, 1989. [12] P.G. Ciarlet.The finite element method for elliptic problems. Amsterdam, North- Holland, 1987. [13] C. Cuvelie, A. Segal and A. van Steenhoven.Finite element methods and Navier- Stokes equations. Reidel, 1986. [14] M. Crouzeix, V. Thomee.The stability in Lpand W1 pof the L2-projection onto finite element function spaces. Math. Comp., 48 , 521-532, 1987. [15] M.Dauge.Elliptic boundary value problems on corner domains. Smoothness and asymptotics of solutions. Lecture Notes in Mathematics, 1341. Springer-Verlag, Berlin, pp. viii+259, 1988. [16] M. Dauge. Stationary Stokes and Navier-Stokes systems on two or three-dimensional domains with corners, SIAM J. Math. Anal. 20 (1), 74-97, 1989. [17] M. Dauge.Probl`emes mixtes pour le laplacien dans des domaines poly´edraux courbes. (French. English summary). C. R. Acad. Sci. Paris S´er. I Math. 309 (8), 553-558, 1989. 50 [18] M. Dauge.Neumann and mixed problems on curvilinear polyhedra, Integr. equat. Oper. Th. 15, 227-261, 1992. [19] J.J. Douglas, T. Dupont, L. Wahlbin. The stability in Lqof the L2-projection into finite element function spaces. Numer. Math. 23, 193-197, 1974/75. [20] P. G. de Gennes. The Physics Of Liquid Crystals. Oxford University Press, London, 1974. [21] P.G. de Gennes, J. Prost. The Physics of Liquid Crystals. Clarendon Press, Oxford, 1995. [22] J. Ericksen. Continuum theory of nematic liquid crystals. Res. Mechanica, 21, 381-392, 1987. [23] J. Ericksen. Conservation laws for liquid crystals. Trans. Soc. Rheol., 5 , 22-34, 1961. [24] A. Ern, J.L. Guermond.Theory and Practice of Finite Elements. Springer–Verlag, 2004. [25] F. C. Frank. On the theory of liquid crystals. Discuss. Faraday Soc. 25, 19-28, 1958. [26] X. Feng, A. Prohl. Analysis of a fully discrete finite element method for the phase field model and approximation of its sharp interface limits. Math. Comp. 73 (246), 541-567, 2004. [27] V. Girault, F. Guill´ en-Gonz´ alez. Mixed formulation, approximation and decoupling algorithm for a nematic liquid crystals model. To appear in Math. Comp. [28] V. Girault, N. Nochetto, R. Scott. Estimates of the finite element Stokes projection in W1,∞.C. R. Math. Acad. Sci. Paris 338 (12), 957-96, 2004. [29] V. Girault, J.L. Lions.Two-grid finite-element schemes for the transient Navier-Stokes problem, M2AN, 35, 945-980, 2001. [30] V. Girault, P.A. Raviart. Finite element methods for Navier-Stokes equations: theory and algorithms. Berlin, Springer-Verlag, 1986. [31] F. Guill´ en-Gonz´ alez, J.V. Guti´ errez-Santacreu. Unconditional stability and convergence of a fully discrete scheme for 2Dviscous fluids models with mass diffusion. Math. Comp. 77 (63), 1495-1524, 2008. [32] F. Guill´ en-Gonz´ alez, J.V.Guti´ errez-Santacreu. Conditional stability and convergence of a fully discrete scheme for 3DNavier-Stokes equations with mass diffusion SIAM J. Numer. Anal. 46 (5), 2276-2308, 2008. 51 [33] F. Guill´ en-Gonz´ alez, J.V. Guti´ errez-Santacreu. A linear mixed finite element scheme for a nematic Ericksen-Leslie liquid crystal model. Submitted. [34] P. Grisvard. Elliptic Problems in Nonsmooth Domains. Pitman Monographs and Studies in Mathematics. 24, Pitman, Boston, MA, 1985. [35] J. G. Heywood, R. Rannacher.Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization. Siam J. Numer. Anal. 19 ( 2), 275–311, 1982. [36] Q. Hu, X.-C. Tai, and R. Winther. A saddle-point approach to the computation of harmonic maps. SIAM J. Numer. Anal. 47 ( 2), 1500-1523, 2009. [37] T.J.R. Hughes, L.P. Franca and M. Balestra.A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuˆska-Brezzi condition: a stable Petrov-Galerkin formulation for the Stokes problem accommodating equal-order interpolations. Computer Methods in Applied Mechanics and Engineering. 59, 85-99, 1986. [38] R. B. Kellogg, J. E. Osborn. A regularity for the Stokes problem in a convex polygon, J. Funct. Anal. 21, 397-431, 1976. [39] F. Leslie. Theory of flow phenomenum in liquid crystals. The Theory of Liquid Crystals. W. Brown, ed. Academic Press, New York, 4, 1-81, 1979. [40] F. Leslie. Some constitutive equations for liquid crystals. Arch. Ration. Mech. Anal., 28 , 265-283, 1968. [41] F.H. Lin. Nonlinear theory of defects in nematic liquid crystals: phase transition and flow phenomena, Comm. Pure Appl. Math. 42, 789-814, 1989. [42] F.H. Lin, C. Liu. Non-parabolic dissipative systems modelling the flow of liquid crystals. Comm. Pure Appl. Math. 48, 501-537, 1995. [43] F. H. Lin, C. Liu. Existence of solutions for the Ericksen-Leslie system. Arch. Rational. Mech. Anal. 154, 135-156, 2000. [44] P. Lin, C. Liu. Simulations of singularity dynamics in liquid crystal flows: A C0finite element approach. Journal of Computational Physics 215, 348-362, 2006. [45] P. Lin, C. Liu, H. Zhang. An energy law preserving C0finite element scheme for simulating the kinematic effects in liquid crystal dynamics. J. Comput. Phys. 227 (2), 1411-1427, 2007. 52 [46] J. L. Lions.On some questions in boundary value problems of mathematics physics. Contempory Developments in Continuum Mechanics and Partial Differential Equations, North- Holland Math. Stud. 30, North-Holland, Amsterdam, 283-346, 1978. [47] J. L. Lions and E. Magenes.Nonhomogeneous Boundary Value Problems and Applications. Vol. I, Translated from the French by P. Kenneth, Die Grundlehren der Mathematischen Wissenschaften, Band 181, Springer-Verlag, New York, 1972. [48] C. Liu, N. J. Walkington.Approximation of liquid crystal flows. SIAM J. Numer. Anal. 37 (3), 725-741, 2000. [49] C. Liu, N. J. Walkington.Mixed methods for the approximation of liquid crystal flows. M2AN Math. Model. Numer. Anal. 36 (2), 205-222, 2002. [50] C. Liu, N. J. Walkington.Convergence of numerical approximations of the incompressible Navier-Stokes equations with variable density and viscosity. SIAM J. Numer. Anal. 45 (3), 1287-1304, 2007. [51] J. Malek, J. Necas, M. Rokyta, M. Ruzicka.Weak and measure-valued solutions to evolutionary PDEs. Chapman & Hall, 1996. [52] J. Necas .Les M´ethodes directes en th´eorie des ´equations elliptiques, Masson, Paris, 1967. [53] C. Oseen. Theory of liquid crystals. Trans. Faraday Soc. 29, 883-899, 1933. [54] M. Protter and H. Weinberger.Maximum Principles in Differential Equations. Prentice-Hall, Englewood Cliffs, NJ, 1967. [55] G. Savar´ e.Regularity results for elliptic equations in Lipschitz domains. J. Funct. Anal. 152 (1), 176-201, 1998. [56] L. R. Scott, S. Zhang. Finite element interpolation of non-smooth functions satisfying boundary conditions. Math. Comp. 54, 483-493, 1990. [57] J. Shatah. Weak solutions and development of singularities in SU(2) σ-model. CPAM 41, 459-469, 1988. [58] J. Simon.Compact sets in the Space Lp(0, T;B). Ann. Mat. Pura Appl., 146, 65-97, 1987. [59] R. Temam. Navier-Stokes equations. North-Holland, 1984. 53 [60] P. A. Raviart, J. M. Thomas. A mixed-finite element method for second order elliptic problems, Mathematical aspects of the finite element method, Lecture Notes in Mathematics. Springer, New York, 1977. 54