Variational and numerical analysis of a -tensor model for smectic-A liquid crystals
Abstract.
We analyse an energy minimisation problem recently proposed for modelling smectic-A liquid crystals. The optimality conditions give a coupled nonlinear system of partial differential equations, with a second-order equation for the tensor-valued nematic order parameter and a fourth-order equation for the scalar-valued smectic density variation . Our two main results are a proof of the existence of solutions to the minimisation problem, and the derivation of a priori error estimates for its discretisation of the decoupled case (i.e., ) using the interior penalty method. More specifically, optimal rates in the and norms are obtained for , while optimal rates in a mesh-dependent norm and norm are obtained for . Numerical experiments confirm the rates of convergence.
Key words and phrases:
interior penalty method, a priori error estimates, finite element methods, smectic liquid crystals1991 Mathematics Subject Classification:
76A15, 49J10, 35B45, 65N12, 65N30Introduction
Smectic liquid crystals (LC) are layered mesophases that have a periodic modulation of the mass density along one spatial direction. Informally, they can be thought of as one-dimensional solids along the direction of periodicity and two-dimensional fluids along the other two remaining directions. According to different in-layer structures, several types of smectic LC are recognised, e.g., smectic-A, smectic-B and smectic-C (see [26, Figure 9] for illustrations of different smectic phases). In particular, in the smectic-A phase, the molecules tend to align parallel to the normals of the smectic layers. For a broader review of liquid crystals, see [2, 16, 28].
Several models have been proposed to describe smectic-A liquid crystals. The classical de Gennes model [15] employs a complex-valued order parameter to describe the magnitude and phase of the density variation. This complex-valued parameterisation leads to some key modelling difficulties, which motivated the development by Pevnyi, Selinger & Sluckin (PSS) of a smectic-A model with a real-valued smectic order parameter and a vector-valued nematic order parameter [24]. The use of a vector-valued nematic order parameter limits the kinds of defects the model can permit [2], and hence Ball & Bedford (BB) proposed a version of the PSS model employing a tensor-valued nematic order parameter [3] instead. The model considered in this work is similar to the BB model, but with additional modifications to render it amenable to numerical simulation (see [31] for details). The model was used to numerically simulate several key smectic defect structures, such as oily streaks and focal conic domains, for the first time.
While the numerical modelling of nematic liquid crystals is now mature, relatively little work has considered smectics. Garcìa-Cervera and Joo [19] perform numerical simulations of the classical de Gennes model [15] in the presence of a magnetic field, using a combined Fourier-finite difference approach. Wittmann et al. use density functional theory to investigate the topology of smectic liquid crystals in complex confinement [30]. Monderkamp et al. examine the topology of defects in two-dimensional confined smectics with the help of extensive Monte Carlo simulations [22]. Our goal in this work is to analyse a model for smectic-A liquid crystals, and its finite element discretization, that was recently proposed by Xia, MacLachlan, Atherton and Farrell in [31].
We consider an open, bounded and convex spatial domain , with Lipschitz boundary . The smectic-A free energy we analyse in this work is given by
| (1) |
where denotes the identity matrix, is the Hessian operator, is the smectic bulk energy density
and is the nematic bulk energy density
| (2) | ||||
Here, is the nematic-smectic coupling parameter, represent smectic bulk constants with , and are nematic elastic and bulk constants, respectively, is the wave number of the smectic layers, and the trace of a matrix is given by .
The two main contributions of this paper are to prove existence of minimisers to the problem of minimising over a suitably-defined admissible set, and to derive a priori error estimates for its discretisation using the interior penalty method [7, 9]. We show the existence result of minimising the free energy eq. 1 in section 1. We then derive a priori error estimates for both and in the simplified case in section 2. We confirm that the theoretical predictions match numerical experiments in section 3, for both and .
In order to avoid the proliferation of constants throughout this work, we use the notation (respectively, ) to represent the relation (respectively, ) for some generic constant (possibly not the same constant on each appearance) independent of the mesh.
1. Existence of minimisers
To formulate the minimisation problem for eq. 1, we must first define the admissible space in which minimisers will be sought. We define as
| (3) |
with the specified Dirichlet boundary data and . Here, denotes the space of symmetric and traceless real-valued matrices. For simplicity of the analysis, we only consider Dirichlet boundary conditions for and here, but other types of boundary conditions (e.g., a mixture of the Dirichlet and natural boundary conditions for as illustrated in the implementations in [31]) can be taken. With this admissible space, we consider the minimisation problem
| (4) |
Notice that is the classical Landau de Gennes (LdG) free energy density [16, 23] for nematic LC and it is proven by Davis & Gartland [14, Corollary 4.4] that there exists a minimiser of the functional for in three dimensions. Moreover, Bedford [4, Theorem 5.18] gives an existence result for the BB model in three dimensions:
with an admissible space , where denotes special functions of bounded variation. For simplicity, we have ignored the surface integral here in the energy functional of the BB model. One can observe its resemblance to eq. 1. Motivated by the above results, we prove the existence of minimisers to eq. 1 via the direct method of calculus of variations.
Theorem 1.1.
Proof.
Note that both the smectic density and the nematic bulk density are bounded from below as . Thus, is also bounded from below and we can choose a minimising sequence , i.e.,
| (5) | ||||
Here we set (resp. ) to be any function with trace (resp. ). We tackle the three terms in separately in the following.
First, for the nematic energy term , we can follow the proof of [14, Theorem 4.3] to obtain that is coercive in , i.e., grows unbounded as , and thus the minimising sequence must be bounded in . Since is a reflexive Banach Space, we have a subsequence (also denoted as ) that weakly converges to such that , and from the Rellich–Kondrachov theorem it follows that The weak lower semi-continuity of the nematic energy density in eq. 2 is guaranteed by [14, Lemma 4.2], therefore,
| (6) |
For the smectic bulk term , we can follow the proof in [4, Theorem 5.19]. By eq. 5, we have
which implies an upper bound for using [4, Ineq. (5.42)]. Hence, is bounded in and thus there is a subsequence (also denoted as ) such that and . Moreover, one can readily check that by the Sobolev embedding of into the Hölder spaces ( for and for ) and the boundedness of domain . Again, it follows from the Rellich–Kondrachov theorem that Noting is bounded from below for all , we then obtain
| (7) |
Now, we consider the nematic-smectic coupling term in . Note that when the wave number , this term reduces to and it is straightforward to obtain the weak lower semi-continuity property. Therefore, we discuss the case of in detail as follows. By the -boundedness property of the minimising sequence and the fact that , we can deduce
Hence, in , and further,
Therefore, we have
| (8) |
Remark 1.2.
We briefly compare the proof with that of Bedford [4, Theorem 5.18]. In that work, the admissible space included a uniaxial constraint.This assumption is necessary to guarantee the -boundedness property of the minimising sequence , since that work seeks instead of . Enforcing this constraint numerically is difficult [6]; in this work we prefer the choice , which enables us to remove the uniaxiality assumption.
2. A priori error estimates
We now consider the discretisation of the minimisation problem eq. 4. For simplicity, we analyse the decoupled case with , where eq. 4 splits into two independent problems: one for the smectic density variation :
and one for the tensor field ,
Here, and . One can derive the following strong forms of their equilibrium equations using integration by parts. The optimality condition for yields a fourth-order partial differential equation (PDE)
where denotes the outward unit normal. Note that the second boundary condition of is a natural boundary condition on the second derivative of . For simplicity, we only consider the case of a cubic nonlinearity (i.e., ) here as the quadratic term can be tackled similarly. Therefore, we consider the following governing equations
| (9) |
Meanwhile, the optimality condition for yields a second-order PDE
| (10) |
We now consider these two problems and in turn.
Remark 2.1.
2.1. A priori error estimates for
Since the PDE eq. 9 for the density variation is a fourth-order problem, a conforming discretisation requires a finite dimensional subspace of the Sobolev space , which necessitates the use of -continuous elements. The construction of these elements is quite involved, particularly in three dimensions; without a special mesh structure, the lowest-degree conforming elements are the Argyris [1] and Zhang [32] elements, of degree 5 and 9 in two and three dimensions respectively. One approach to avoid such complexity is to use mixed formulations by solving two second-order systems, and we refer to [12, 27] for instance. However, this substantially increases the size of the linear systems to be solved. Alternatively, one can directly tackle the fourth-order problem with non-conforming elements, that do not satisfy the -requirement. For instance, the so-called continuous/discontinuous Galerkin methods and interior penalty methods (-IP) are analysed in [9, 17], combining concepts from the theory of continuous and discontinuous Galerkin methods. Essentially, these methods use -conforming elements and penalise inter-element jumps in first derivatives to weakly enforce -continuity. This has the advantages of both convenience and efficiency: the weak form is simple, with only minor modifications from a conforming method, and fewer degrees of freedom are used than with a fully discontinuous Galerkin method.
We thus adopt the idea of -IP methods to solve the nonlinear fourth-order problem and derive a priori error estimates regarding . We adapt the techniques of [21] to prove convergence rates with the use of familiar continuous Lagrange elements for the problem . The weak form of eq. 9 is defined as: find such that
| (11) |
where for ,
and for ,
Since eq. 11 is nonlinear, we derive its linearisation: find such that
| (12) |
where represents the dual pairing between and .
It is straightforward to derive the coercivity and boundedness of the bilinear operator with the semi-norm (in fact, this is indeed a norm in ).
Lemma 2.2.
For , there holds
Let be a mesh of with denoting an element, and let (resp. ) denote the set of all interior (resp. boundary) edges/faces of the mesh, and . Define the broken Sobolev space equipped with the broken norm . We take a nonconforming but continuous approximation for the solution of eq. 11, that is to say, defined for (since is a fourth-order problem) by
with denoting piecewise continuous polynomials of degree on a mesh of quadrilaterals or hexahedra. Following the derivation of the -IP formulation given in [7, Section 3], we introduce the discrete nonlinear weak form: find such that
| (13) |
where for all ,
and
| (14) |
Here, is the penalty parameter (to be specified in section 3), indicates the size of the edge/face and the average of the second derivatives of along the normal direction across the edge/facet is defined as For any interior edge shared by two cells and , we define the jump by with representing the restriction of outward normals in respectively. On the boundary edge/face , we define . The operator penalises the first derivatives across the interior edge/facet since the function in is not necessarily continuously differentiable.
The linearised version is to seek such that
| (15) |
where
| (16) |
We also define the mesh-dependent -like semi-norm for ,
| (17) |
Note that is indeed a norm on . This norm will be used in the well-posedness and convergence analysis below.
We first give an immediate result about the consistency of the discrete form eq. 13.
Theorem 2.3.
Proof.
Multiplying the fourth-order term in eq. 9 with a test function and using piecewise integration by parts with the boundary condition specified in eq. 9 for , one can obtain
| (18) |
Since implies is continuous on the whole domain , the jump term then becomes zero and we can thus symmetrise and penalise the form eq. 18. This leads to the presence of . The remaining terms involving and are straightforward as one takes the test function . Therefore, satisfies eq. 13. ∎
By noting that it is natural to extend the classical elliptic regularity result [5] for the biharmonic operator to the case of the bi-Hessian operator . If the domain is smooth, the weak solutions of the biharmonic problem (e.g., [7, Example 2]) belong to by classical elliptic regularity results and thus we make this assumption henceforth.
Moreover, to facilitate the analysis, we further assume that is an isolated solution, i.e., there is only one solution satisfying eq. 9 within a sufficiently small ball with radius . These assumptions then imply that the linearised operator satisfies the following inf-sup condition:
| (19) |
2.1.1. Well-posedness of the discrete form
Recalling [7, Eq. (3.20)], we can obtain for ,
| (20) |
as the edge/facet size . With the estimate eq. 20 at hand, we then apply the Cauchy–Schwarz inequality and use the definition eq. 17 of to obtain the boundedness of and . We omit the details of the proofs here and only illustrate the boundedness result for and below.
Lemma 2.4.
(Boundedness of and .) For , we have
| (21) |
For , ,
| (22) |
Proof.
By Hölder’s inequality, the Sobolev embedding , and the fact that the semi-norm is a norm in , we deduce
It then follows from a Poincaré inequality [11, Eq. (5.7)] for piecewise functions that
| (23) |
Thus, we obtain
We give the coercivity result for the bilinear form .
Lemma 2.5.
(Coercivity of ) For a sufficiently large penalty parameter , there holds
| (24) |
Proof.
An important question about the well-posedness is the coercivity of the bilinear operator . Due to the presence of and terms in , it is not trivial to derive its coercivity. We first discuss the weak coercivity of the bilinear form defined as
| (25) |
To this end, we will employ the enrichment operator with being the Hsieh–Clough–Tocher macro finite element space [7]. The following lemma is adapted to our notation and definition of from [7, Lemma 1].
Lemma 2.6.
[7, Lemma 1] For , there holds
With this, we obtain the following discrete inf-sup condition for the discrete bilinear operator .
Theorem 2.7.
(Weak coercivity of ) Let be a regular isolated solution of the nonlinear continuous weak form eq. 13. For a sufficiently large and a sufficiently small mesh size , the following discrete inf-sup condition holds with a positive constant :
| (26) |
Proof.
For , it follows from the boundedness result of that and . Furthermore, since is bounded and coercive as given by lemma 2.2, for a given with , there exists and that solve the linear systems:
| (27a) | |||
| (27b) |
We require the -regularity for the derivation of eq. 34 below. It then follows from the standard elliptic regularity result that with a constant depending on .
Subtracting eq. 27a from eq. 27b, then taking and using the coercivity of , boundedness of , and enrichment estimate lemma 2.6, we get
| (28) |
Since is a regular isolated solution of eq. 11, it yields by eq. 19, eq. 27b, lemma 2.2, the fact that and the triangle inequality, that there exists with such that
| (29) |
Note that on since . We can thus calculate using lemma 2.6 that
Further, by the triangle inequality, we get
| (30) |
Since , it follows from the coercivity result eq. 24 that there exists with such that
| (31) |
where in the last equality we have used the fact that
because of eq. 27a and .
Using the boundedness result lemma 2.4 and the enrichment estimate lemma 2.6, we obtain
| (32) |
By the boundedness of the bilinear form and standard interpolation estimate, we have
| (33) |
Moreover, by the fact that , the enrichment estimate lemma 2.6 and eq. 18, there holds
| (34) |
Combine eqs. 32, 33 and 34 in eq. 31 to obtain
| (35) |
Substituting eq. 35 into eq. 30 and using standard interpolation estimates yield that
| (36) |
A use of eqs. 36 and 35, standard interpolation estimates and eq. 28 in eq. 29 leads to
Then, by the triangle inequality, we have
Therefore, for the mesh size satisfying the discrete inf-sup condition eq. 26 holds for . ∎
We can now obtain the discrete inf-sup condition for the perturbed bilinear form given by
| (37) |
Here, we employ the interpolation operator .
Theorem 2.8.
(Weak coercivity of ) Let be a regular isolated solution of the nonlinear continuous weak form eq. 13 and the interpolation of . For a sufficiently large and a sufficiently small mesh size , the following discrete inf-sup condition holds:
| (38) |
Proof.
Denote . By the definition eq. 37 of the bilinear form , we have It follows from the definition of and its boundedness result lemma 2.4 that
where is the generic constant arising in the boundedness result lemma 2.4 for . Therefore, we obtain that
Now using the inf-sup condition theorem 2.7 for the bilinear form , boundedness result lemma 2.4 and interpolation estimates, we get
for a sufficiently small mesh size such that . Here, depends on and and gives the regularity of , i.e., . Therefore, the inf-sup condition eq. 38 holds. ∎
2.1.2. Convergence analysis
We proceed to the error analysis for the discrete nonlinear problem eq. 13. Let We define the nonlinear map by
| (39) |
for . Due to the weak coercivity property in theorem 2.8, the nonlinear map is well-defined.
The existence and local uniqueness of the solution to the discrete nonlinear problem eq. 13 will be proven via an application of Brouwer’s fixed point theorem, which necessitates the use of two auxiliary lemmas illustrating that (i) maps from a ball to itself; and (ii) the map is contracting.
Lemma 2.9.
(Mapping from a ball to itself) Let be a regular isolated solution of the continuous nonlinear weak problem eq. 11. For a sufficiently large and a sufficiently small mesh size , there exists a positive constant such that:
Proof.
Note that the solution of eq. 11 satisfies the discrete weak formulation eq. 13 due to the consistency result theorem 2.3, that is to say, there holds that
| (40) |
By the linearity of , the definition eq. 39 of the nonlinear map and eq. 40, we calculate
In what follows, we give upper bounds for each . Using the boundedness of and the interpolation estimate [9, Eq. (5.3)] in the -norm, we obtain
We rearrange terms in and use the boundedness result lemma 2.4 and the interpolation result [9, Eq. (5.3)] to obtain
Let . We use the definition of and its boundedness to deduce that
Hence, we combine the above bounds for , to have
By the inf-sup condition eq. 38 for the perturbed bilinear form, we further deduce that there exists a with such that Since , we obtain
Note that there are other cases when and we only focus on the case of here for brevity. The idea of the remainder of the proof is to choose an appropriate so that . For simplicity of the calculation, we illustrate the case when . To this end, we take and choose satisfying This yields
This completes the proof. ∎
Lemma 2.10.
(Contraction result) For a sufficiently large , a sufficiently small mesh size and any , there holds
| (41) |
Proof.
For , we use the definition eq. 39 of the nonlinear map , the definition eq. 37 and linearity of to calculate
Let , and . We make some elementary manipulations and use the boundedness of and the inequality of geometric and arithmetic means to yield
By the inf-sup condition eq. 38, we know that there exists with such that
Therefore, we have Noting that for completes the proof. ∎
The existence and local uniqueness of the discrete solution can now be obtained via the application of Brouwer’s fixed point theorem [20].
Theorem 2.11.
Proof.
A use of lemma 2.9 yields that the nonlinear map maps a closed convex set to itself. Moreover it is a contracting map. Therefore, an application of Brouwer’s fixed point theorem yields that has at least one fixed point, say , in this ball . The uniqueness of the solution to eq. 13 in that ball follows from the contraction result in lemma 2.10. Meanwhile, we have by lemma 2.9 that
| (42) |
The error estimate is then obtained straightforwardly using the triangle inequality combined with eq. 42 and the interpolation estimate [9, Eq. (5.3)]. ∎
It follows from theorem 2.11 that optimal convergence rates are achieved in the mesh-dependent norm . This will be numerically verified in section 3.
2.1.3. Estimates in the -norm
We derive an error estimate using a duality argument in this subsection. To this end, we consider the following linear dual problem to the primal nonlinear problem eq. 9:
| (43) |
for . For smooth domains , it can be deduced by a classical elliptic regularity result that . The corresponding weak form is: find such that
that is to say,
| (44) |
Remark 2.12.
The first equality in eq. 44 holds since and .
We give two auxiliary results in the following.
Lemma 2.13.
For , , and , there holds that
Proof.
Note that since and on . We calculate
Here, the last, second last, and third last steps follow from standard interpolation estimates, the Cauchy–Schwarz inequality, and integration by parts twice, respectively. ∎
Lemma 2.14.
The solution of the linear dual problem eq. 43 belongs to on a smooth domain and it holds that
| (45) |
Proof.
We can use the inf-sup condition eq. 19 for the linear operator , the weak form eq. 44 and the Cauchy–Schwarz inequality to obtain
| (46) |
By the form of eq. 44, the boundedness of and , and eq. 46, we have
| (47) |
Using a bootstrapping argument in elliptic regularity (see, e.g., [18, Section 6.3]), we can deduce that in a smooth domain . The regularity estimate eq. 45 follows from eq. 47. ∎
We are ready to derive the a priori error estimates.
Theorem 2.15.
( error estimate) Under the same conditions as theorem 2.11 and assuming further that (since the problem is fourth-order), the discrete solution approximates such that
| (48) |
Proof.
Taking in eq. 43 and multiplying eq. 43 by a test function and integrating by parts, we obtain It follows from the fact that , , and the definition eq. 11 of the nonlinear continuous weak form that
We bound each separately using the boundedness of and , theorem 2.11 and standard interpolation estimates. This leads to
Setting and estimating as in of lemma 2.9 with the use of theorem 2.11 and yields
Combining the above estimates for () and using the regularity estimate eq. 45 and , we obtain
Hence, eq. 48 follows from the triangle inequality and standard interpolation estimates. ∎
theorem 2.15 implies that for quadratic approximations to the sufficiently regular solution of eq. 9, there is a sub-optimal convergence rate in the -norm while for higher order () approximations, we expect optimal error rates. We shall see numerical verifications of this in the subsequent sections.
2.1.4. The inconsistent discrete form
The above analysis considers the consistent weak formulation eq. 13. In practice, Xia et al. [31] adopted the inconsistent discrete weak form in the implementation due to its cheaper assembly cost: find such that
| (49) |
where Comparing and , the missing terms are the interior facet integrals arising from piecewise integration by parts and symmetrisation. Due to the absence of these terms in , one can immediately notice that the discrete weak formulation eq. 49 is inconsistent in the sense that the solution of the strong form eq. 9 does not satisfy the weak form eq. 49, as opposed to the result of theorem 2.3.
Despite this inconsistency, in practice this also leads to a convergent numerical scheme with similar convergence rates, as illustrated in section 3. This is not surprising; a similar idea has also been applied and introduced as weakly over-penalised symmetric interior penalty (WOPSIP) methods in [10] for second-order elliptic PDEs and in [8] for biharmonic equations.
Remark 2.16.
The excessive size of the penalty parameter in the WOPSIP method could induce ill-conditioned linear systems. No such effects are observed in our numerical results.
2.2. A priori error estimates for
Problem is a special form of the classical LdG model of nematic LC. Finite element analysis for a more general form using conforming discretisations has been studied in [13, 14]. More specifically, Davis and Gartland [14] gave an abstract nonlinear finite element convergence analysis where an optimal error bound is proved on convex domains with piecewise linear polynomial approximations, but do not derive an error bound in the norm. Recently, Maity, Majumdar and Nataraj [21] analysed the discontinuous Galerkin finite element method for a two-dimensional reduced LdG free energy, where optimal a priori error estimates in the -norm with exact solutions in are achieved for a piecewise linear discretisation. Both works only focus on piecewise linear approximations. In this section, we will follow similar steps to section 2.1 to prove the - and -convergence rates for the problem with the use of common continuous Lagrange elements of arbitrary positive degree. Since the approach is similar to the previous subsections, we omit some details for brevity.
The continuous weak formulation of in two dimensions (the three-dimensional case can be tackled similarly) is given by: find such that
| (50) |
where the bilinear forms are and the nonlinear operator is given by
| (51) |
Since eq. 50 is nonlinear, we need to approximate the solution of its linearised version, i.e., find such that
| (52) |
where represents the dual pairing between and .
Suppose approximates the solution of eq. 50 with the conforming finite element method on a finite dimensional space . Throughout this subsection we take . Furthermore, we denote and . We assume that the minimiser to be approximated is isolated, i.e., the linearised operator is nonsingular. This is equivalent to the following continuous inf-sup condition [21, Eq. (2.8)]:
| (53) |
With this inf-sup condition for , we can obtain a stability result for the perturbed bilinear form by following similar steps as in the proof of theorem 2.8.
Theorem 2.17.
(Stability of the perturbed bilinear form) Let be a regular isolated solution of the nonlinear continuous weak form eq. 50 and the interpolant of . For a sufficiently small mesh size , the following discrete inf-sup condition holds:
| (54) |
We give some auxiliary results about the operators , and that can be verified via the Cauchy–Schwarz inequality, the Poincaré inequality, and Sobolev embeddings.
Lemma 2.18.
(Boundedness and coercivity of ) For , there holds
Lemma 2.19.
(Boundedness of , ) For , there holds
| (55) |
and for , , ,
| (56) |
To proceed to error estimates for the nonlinear problem eq. 50, we define the nonlinear map by
for . Due to the stability result of theorem 2.17, the nonlinear map is well-defined. We define the local ball . The following two auxiliary lemmas provide the necessary components for the application of Brouwer’s fixed point theorem.
Lemma 2.20.
(Mapping from a ball to itself) Let be a regular isolated solution of the continuous nonlinear weak problem eq. 50. For a sufficiently small mesh size , there exists a positive constant such that:
Remark 2.21.
In fact, the choice of can be taken as in the proof of lemma 2.20. Here, denotes the regularity index of , i.e., .
Lemma 2.22.
(Contraction result) For a sufficiently small mesh size and any , there holds
| (57) |
Remark 2.23.
In the proof of lemma 2.22, we have particularly used the stability property of the perturbed bilinear form as given by theorem 2.17.
Hence, the existence and local uniqueness of the discrete solution can be derived by following similar steps as in the proof of theorem 2.11.
Theorem 2.24.
We again employ an Aubin–Nitsche duality argument to derive error estimates. To this end, we consider the following linear dual problem to the primal nonlinear problem eq. 10: find such that
| (58) |
for a given (we will make a particular choice for in the proof of theorem 2.27). The weak form of eq. 58 is to find such that
| (59) |
To derive the a priori error estimates, we need two more auxiliary results.
Lemma 2.25.
For , , and , it holds that
Lemma 2.26.
(Boundedness of the dual solution in the -norm) The solution to the weak form eq. 59 of the dual linear problem belongs to and it holds that
| (60) |
Finally, we are ready to deduce an optimal error estimate.
Theorem 2.27.
We will verify these results in the next section.
3. Numerical experiments
The proceeding section presents some a priori error estimates for both and in the decoupled case . We now test the convergence rate of the finite element approximations by the method of manufactured solutions (MMS) and experimentally investigate the coupled case in two dimensions. To this end, we choose a nontrivial solution for each state variable and add an appropriate source term to the equilibrium equations, thus modifying the energy accordingly. We can then compute the numerical convergence order.
3.1. Test 1: on the unit square
In this test, the numerical runs are performed on the unit square and we take the following exact expressions for each state variable,
| (62) | ||||
Then, in conducting the MMS, we are to solve the following governing equations
subject to Dirichlet boundary conditions for both and and a natural boundary condition for . Here, source terms , and are derived by substituting eq. 62 to the left hand sides, and and are given by
We partition the domain into squares with uniform mesh size (, , , ) and denote the numerical solutions by , and . The numerical errors of and in the -, - and -norms are defined as
The convergence order is then calculated from the formula Throughout this section, we use the parameter values , , , , and , similar to the simulations of oily streaks in [31].
Remark 3.1.
Since this is purely a numerical verification exercise, the manufactured solution can be physically unrealistic. However, we must specify a reasonable initial guess for Newton’s method, due to the nonlinearity of the problem. The initial guess throughout this section is taken to be .
3.1.1. Convergence rate for
For we expect both optimal and rates, as illustrated in theorems 2.27 and 2.24. table 1 presents the numerical convergence rate for the finite elements , and . Optimal and rates are shown with all choices of finite elements, as predicted.
| rate | rate | ||||
| 6 | 8.12 | – | 3.78 | – | |
| 12 | 2.02 | 2.01 | 1.88 | 1.01 | |
| 24 | 5.05 | 2.00 | 9.39 | 1.00 | |
| 48 | 1.26 | 2.00 | 4.69 | 1.00 | |
| 6 | 2.92 | – | 1.11 | – | |
| 12 | 3.90 | 2.90 | 2.71 | 2.04 | |
| 24 | 5.02 | 2.96 | 6.72 | 2.01 | |
| 48 | 6.36 | 2.99 | 1.68 | 2.00 | |
| 6 | 3.02 | – | 2.25 | – | |
| 12 | 2.17 | 3.80 | 2.72 | 3.05 | |
| 24 | 1.45 | 3.90 | 3.34 | 3.03 | |
| 48 | 9.33 | 3.96 | 4.13 | 3.01 |
Regarding the density variation , we first present the convergence behaviour of the consistent discrete formulation eq. 13 with penalty parameter , since we have proven the optimal error rate in the mesh-dependent norm . The errors and convergence orders are listed in table 2. Optimal rates are observed in the -norm. Furthermore, optimal orders of convergence in the -norm are shown for approximating polynomials of degree greater than , while a sub-optimal rate in the -norm is given for piecewise quadratic polynomials, exactly as expected. Sub-optimal convergence rates for quadratic polynomials were also illustrated in the numerical results of [29]. We also tested the convergence with the penalty parameter and found that the discrete norms are very similar to table 2. We therefore avoid repeating the details here.
| rate | rate | rate | |||||
|---|---|---|---|---|---|---|---|
| 6 | 1.17 | – | 3.46 | – | 1.36 | – | |
| 12 | 2.60 | 2.17 | 9.81 | 1.82 | 7.25 | 0.91 | |
| 24 | 6.37 | 2.03 | 2.54 | 1.95 | 3.54 | 1.03 | |
| 48 | 1.82 | 1.80 | 6.88 | 1.88 | 1.76 | 1.01 | |
| 6 | 4.73 | – | 1.32 | – | 4.98 | – | |
| 12 | 3.32 | 3.83 | 1.41 | 3.23 | 9.96 | 2.32 | |
| 24 | 2.12 | 3.97 | 1.63 | 3.12 | 2.46 | 2.02 | |
| 48 | 1.32 | 4.00 | 1.99 | 3.03 | 6.14 | 2.00 | |
| 6 | 2.01 | – | 7.76 | – | 3.94 | – | |
| 12 | 5.40 | 5.22 | 4.30 | 4.17 | 4.88 | 3.01 | |
| 24 | 1.68 | 5.00 | 2.68 | 4.00 | 6.11 | 2.99 | |
| 48 | 5.27 | 4.99 | 1.68 | 3.99 | 7.64 | 3.00 |
We next give the error rates for the inconsistent discrete formulation eq. 49. We illustrate the discrete norms and the computed convergence rates in table 3 with the penalty parameter . It can be observed that only first order convergence is obtained in the -like norm even with different approximating polynomials. Moreover, we notice by comparing tables 2 and 3 that the convergence rate deteriorates slightly for polynomials of degree 3 (although not for degree 4). This, however, can be improved by choosing a larger penalty parameter, as shown in table 4 with , where optimal rates are shown for the discrete norms , and for all polynomial degrees (except only sub-optimal in when a piecewise quadratic polynomial is used as the approximation). The inconsistent discrete formulation appears to be a reasonable choice when a sufficiently large penalty parameter is used.
| rate | rate | rate | |||||
|---|---|---|---|---|---|---|---|
| 6 | 3.50 | – | 1.06 | – | 5.60 | – | |
| 12 | 8.76 | 5.32 | 5.41 | 4.29 | 2.56 | 1.13 | |
| 24 | 1.77 | 2.31 | 7.47 | 2.86 | 1.28 | 0.99 | |
| 48 | 4.35 | 2.02 | 1.24 | 2.56 | 6.42 | 1.00 | |
| 6 | 6.47 | – | 1.86 | – | 7.59 | – | |
| 12 | 3.40 | 4.25 | 1.73 | 3.43 | 2.74 | 1.47 | |
| 24 | 1.98 | 4.10 | 2.03 | 3.09 | 1.31 | 1.07 | |
| 48 | 3.73 | 2.39 | 2.63 | 2.95 | 6.45 | 1.02 | |
| 6 | 2.05 | – | 7.85 | – | 3.93 | – | |
| 12 | 5.40 | 5.24 | 4.31 | 4.19 | 4.88 | 3.01 | |
| 24 | 1.68 | 5.00 | 2.68 | 4.01 | 6.11 | 3.00 | |
| 48 | 5.27 | 5.00 | 1.67 | 4.00 | 7.64 | 3.00 |
| rate | rate | rate | |||||
|---|---|---|---|---|---|---|---|
| 6 | 1.17 | – | 3.48 | – | 1.36 | – | |
| 12 | 2.62 | 2.16 | 9.86 | 1.82 | 7.26 | 0.91 | |
| 24 | 6.38 | 2.04 | 2.54 | 1.96 | 3.54 | 1.03 | |
| 48 | 1.82 | 1.81 | 6.88 | 1.88 | 1.76 | 1.01 | |
| 6 | 4.80 | – | 1.35 | – | 4.92 | – | |
| 12 | 3.35 | 3.84 | 1.43 | 3.23 | 9.86 | 2.32 | |
| 24 | 2.14 | 3.97 | 1.63 | 3.13 | 2.45 | 2.01 | |
| 48 | 1.33 | 4.01 | 1.99 | 3.04 | 6.13 | 2.00 | |
| 6 | 2.05 | – | 7.85 | – | 3.93 | – | |
| 12 | 5.40 | 5.24 | 4.31 | 4.19 | 4.88 | 3.01 | |
| 24 | 1.68 | 5.00 | 2.68 | 4.01 | 6.11 | 3.00 | |
| 48 | 5.27 | 5.00 | 1.67 | 4.00 | 7.64 | 3.00 |
3.1.2. Convergence rate for
We next investigate the numerical convergence behaviour in the coupled case, i.e., , in this subsection. Its analysis remains future work, but since it is the coupled case that is solved in practice it is important to assure ourselves that the discretisation is sensible. For brevity, we fix the model parameter .
We examine the inconsistent discretisation for with penalty parameter . In unreported preliminary experiments, we observed that the error in is governed by the lower of the degrees of the polynomials used for and . We thus give the convergence rates for and separately in tables 5 and 6, with the other degree fixed appropriately. It can be seen that retains optimal rates in both the and norms, and though there are some fluctuations of the order for , it still possesses very similar convergence rates when compared with the decoupled case described in table 4.
| rate | rate | rate | |||||
|---|---|---|---|---|---|---|---|
| 6 | 1.21 | – | 3.59 | – | 1.37 | – | |
| 12 | 3.98 | 1.61 | 1.42 | 1.34 | 8.30 | 0.72 | |
| 24 | 1.57 | 1.35 | 4.99 | 1.51 | 3.89 | 1.09 | |
| 48 | 2.58 | 2.60 | 9.06 | 2.46 | 1.78 | 1.13 | |
| 6 | 7.36 | – | 2.25 | – | 9.10 | – | |
| 12 | 4.13 | 4.16 | 1.86 | 3.60 | 1.11 | 3.03 | |
| 24 | 4.23 | 3.29 | 2.24 | 3.05 | 2.53 | 2.14 | |
| 48 | 3.01 | 3.81 | 2.28 | 3.29 | 6.15 | 2.04 |
| rate | rate | ||||
|---|---|---|---|---|---|
| 6 | 8.12 | – | 3.78 | – | |
| 12 | 2.02 | 2.01 | 1.88 | 1.01 | |
| 24 | 5.05 | 2.00 | 9.39 | 1.00 | |
| 48 | 1.26 | 2.00 | 4.69 | 1.00 | |
| 6 | 2.92 | – | 1.11 | – | |
| 12 | 3.90 | 2.90 | 2.71 | 2.04 | |
| 24 | 5.02 | 2.96 | 6.72 | 2.01 | |
| 48 | 6.37 | 2.98 | 1.68 | 2.00 | |
| 6 | 3.02 | – | 2.25 | – | |
| 12 | 2.17 | 3.80 | 2.72 | 3.05 | |
| 24 | 1.45 | 3.90 | 3.34 | 3.03 | |
| 48 | 9.32 | 3.96 | 4.13 | 3.01 |
3.2. Test 2: on the unit disc
For the second set of experiments, we provide numerical results for an exact solution with only -regularity, instead of as in eq. 62. Our goal is to investigate whether the -regularity assumption on can be relaxed. To this end, we consider a triangular mesh of and choose to be
| (63) |
and choose the same exact solution for and as in eq. 62. The exact solution given by eq. 63 is in but not in , hence violating the regularity assumption of the analysis in section 2.1.
The resulting convergence rates are reported in tables 7 and 8. Table 7 shows that optimal and rates are achieved for with three different choices of finite elements . Table 8 shows the convergence behaviour with penalty parameter when using the inconsistent discrete formulation eq. 49. In contrast to table 3, only first order convergence is obtained for the discrete norm and second-order convergence for and , with both and . Interestingly, table 8 indicates no convergence when using elements. It appears that the assumption is necessary for our analysis, and that a different analysis should be carried out when this assumption no longer holds.
| No. of triangles | rate | rate | |||
|---|---|---|---|---|---|
| 60 | 6.08 | – | 1.09 | – | |
| 240 | 1.56 | 1.96 | 5.80 | 0.91 | |
| 960 | 3.92 | 2.00 | 2.93 | 0.99 | |
| 3840 | 9.83 | 2.00 | 1.47 | 1.00 | |
| 15360 | 2.47 | 1.99 | 7.34 | 1.00 | |
| 60 | 8.97 | – | 2.11 | – | |
| 240 | 1.51 | 2.57 | 5.87 | 1.84 | |
| 960 | 2.22 | 2.77 | 1.52 | 1.95 | |
| 3840 | 3.02 | 2.88 | 3.85 | 1.98 | |
| 15360 | 3.93 | 2.94 | 9.67 | 1.99 | |
| 60 | 1.08 | – | 3.21 | – | |
| 240 | 8.21 | 3.72 | 4.58 | 2.81 | |
| 960 | 5.52 | 3.89 | 5.92 | 2.95 | |
| 3840 | 3.54 | 3.96 | 7.44 | 2.99 | |
| 15360 | 2.23 | 3.99 | 9.31 | 3.00 |
| No. of triangles | rate | rate | rate | ||||
|---|---|---|---|---|---|---|---|
| 60 | 1.36 | – | 2.77 | – | 7.39 | – | |
| 240 | 1.58 | 3.11 | 3.86 | 2.84 | 2.39 | 1.63 | |
| 960 | 7.76 | 1.02 | 1.57 | 1.30 | 1.44 | 0.73 | |
| 3840 | 1.84 | -1.25 | 3.42 | -1.12 | 1.26 | 0.19 | |
| 15360 | 2.77 | -0.59 | 5.29 | -0.63 | 1.60 | -0.34 | |
| 60 | 9.94 | – | 2.21 | – | 6.15 | – | |
| 240 | 3.99 | 1.32 | 8.49 | 1.38 | 3.03 | 1.02 | |
| 960 | 1.33 | 4.90 | 4.57 | 4.22 | 1.27 | 1.26 | |
| 3840 | 3.21 | 2.06 | 8.47 | 2.43 | 0.66 | 0.93 | |
| 15360 | 9.22 | 1.80 | 2.11 | 2.00 | 0.34 | 0.97 | |
| 60 | 7.17 | – | 1.65 | – | 4.70 | – | |
| 240 | 1.34 | 2.42 | 3.84 | 2.10 | 2.39 | 0.98 | |
| 960 | 1.18 | 3.50 | 4.54 | 3.08 | 1.28 | 0.90 | |
| 3840 | 2.98 | 1.99 | 8.14 | 2.48 | 0.67 | 0.94 | |
| 15360 | 8.15 | 1.87 | 1.86 | 2.13 | 0.34 | 0.97 |
References
- [1] J. H. Argyris, I. Fried, and D. W. Scharpf, The TUBA family of plate elements for the matrix displacement method, Aeronaut. J., 72 (1968), pp. 701–709.
- [2] J. M. Ball, Mathematics and liquid crystals, Mol. Cryst. Liq. Cryst., 647 (2017), pp. 1–27.
- [3] J. M. Ball and S. J. Bedford, Discontinuous order parameters in liquid crystal theories, Mol. Cryst. Liq. Cryst., 612 (2015), pp. 1–23.
- [4] S. J. Bedford, Calculus of variations and its application to liquid crystals, PhD thesis, University of Oxford, 2014.
- [5] H. Blum and R. R. Bonn, On the boundary value problem of the biharmonic operator on domains with angular corners, Math. Mech. in the Appli. Sci., 2 (1980), pp. 556–581.
- [6] J. P. Borthagaray, R. H. Nochetto, and S. W. Walker, A structure-preserving FEM for the uniaxially constrained Q-tensor model of nematic liquid crystals, Numer. Math., 145 (2020), pp. 837–881.
- [7] S. C. Brenner, interior penalty methods, in Frontiers in Numerical Analysis - Durham 2010. Lecture Notes in Computational Science and Engineering, J. Blowey and M. Jensen, eds., vol. 85, Springer, Berlin, Heidelberg, 2011.
- [8] S. C. Brenner, T. Gudi, and L. Sung, A weakly over-penalized symmetric interior penalty method for the biharmonic problem, Electon. Trans. Numer. Anal., 37 (2010), pp. 214–238.
- [9] S. C. Brenner and L. Sung, interior penalty methods for fourth order elliptic boundary value problems on polygonal domains, J. Sci. Comput., 22 (2005), pp. 83–118.
- [10] S. C. Brenner and L. Sung, A weakly over-penalized symmetric interior penalty method, Electon. Trans. Numer. Anal., 30 (2008), pp. 107–127.
- [11] S. C. Brenner, K. Wang, and J. Zhao, Poincaré-Friedrichs inequalities for piecewise functions, Numer. Funct. Anal. Optim., 25 (2004), pp. 463–478.
- [12] X. Cheng, W. Han, and H. Huang, Some mixed finite element methods for the biharmonic equation, J. Comp. Appl. Math., 126 (2000), pp. 91–109.
- [13] T. A. Davis, Finite element analysis of the Landau–de Gennes minimization problem for liquid crystals in confinement, PhD thesis, Kent State University, 1994.
- [14] T. A. Davis and J. E. C. Gartland, Finite element analysis of the Landau-de Gennes minimization problem for liquid crystals, SIAM J. Numer. Anal., 35 (1998), pp. 336–362.
- [15] P. G. de Gennes, An analogy between superconductors and smectic A, Solid State Commun., 10 (1972), pp. 753–756.
- [16] P. G. de Gennes, The Physics of Liquid Crystals, Oxford University Press, Oxford, 1974.
- [17] G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, and R. L. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 3669–3750.
- [18] L. C. Evans, Partial Differential Equations, vol. 19 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, 2nd ed., 2010.
- [19] C. J. García-Cervera and S. Joo, Layer undulations in Smectic A liquid crystals, J. Comp. Theo. Nano., 7 (2010), pp. 795–801.
- [20] S. Kesavan, Topics in Functional Analysis and Applications, John Wiley & Sons, New York, 1989.
- [21] R. R. Maity, A. Majumdar, and N. Nataraj, Discontinuous Galerkin finite element methods for the Landau–de Gennes minimization problem of liquid crystals, IMA J. Numer. Anal., 41 (2021), pp. 1130–1163.
- [22] P. Monderkamp, R. Wittmann, L. B. G. Cortes, D. G. A. L. Aarts, and H. Löwen, Topology of orientational defects in confined smectic liquid crystals, Phys. Rev. Lett., 127 (2021), pp. 1–10.
- [23] N. J. Mottram and C. J. P. Newton, Introduction to -tensor theory, 2014, https://arxiv.org/abs/1409.3542.
- [24] M. Y. Pevnyi, J. Selinger, and T. J. Sluckin, Modeling smectic layers in confined geometries: order parameter and defects, Phys. Rev. E, 90 (2014), pp. 1–8.
- [25] M. Robinson, C. Luo, P. E. Farrell, R. Erban, and A. Majumdar, From molecular to continuum modelling of bistable liquid crystal devices, Liq. Cryst., 44 (2017), pp. 2267–2284.
- [26] H. Sackmann, Plenary lecture. Smectic liquid crystals. A historical review, Liq. Cryst., 5 (1989), pp. 43–55.
- [27] R. Scholtz, A mixed method for fourth-order problems using linear finite elements, RAIRO Numer. Anal., 15 (1978), pp. 85–90.
- [28] I. W. Stewart, The Static and Dynamic Continuum Theory of Liquid Crystals: A Mathematical Introduction, CPC Press, 2004.
- [29] E. Süli and I. Mozolevski, -version interior penalty DGFEMs for the biharmonic equation, Comput. Methods Appl. Mech. Engrg., 196 (2007), pp. 1851–1863.
- [30] R. Wittmann, L. B. G. Cortes, H. Löwen, and D. G. A. L. Aarts, Particle-resolved topological defects of smectic colloidal liquid crystals in extreme confinement, Nat. Comm., 12 (2021), pp. 1–10.
- [31] J. Xia, S. MacLachlan, T. J. Atherton, and P. E. Farrell, Structural landscapes in geometrically frustrated smectics, Phys. Rev. Lett., 126 (2021), pp. 1–6.
- [32] S. Zhang, A family of 3D continuously differentiable finite elements on tetrahedral grids, Appl. Numer. Math., 59 (2009), pp. 219–233.