Convergence of a fully discrete finite element method for the Beris-Edwards system of liquid crystal dynamics

1. Introduction

• The Landau–de Gennes Q-tensor is a symmetric, traceless second-order tensor used to describe defects that vector-based models cannot represent in all physically occurring forms.The Beris–Edwards model couples incompressible Navier–Stokes equations to Q-tensor gradient flow.Source: §1 · p. 1
• As h and Δt tend to zero with h^2 = o(Δt), the scheme converges along a subsequence to a weak solution; this condition is needed for convergence, not solvability or stability.Weak-solution uniqueness is open, so the result is subsequential.Source: §1 · pp. 3–4
• The fully discrete method combines formally second-order BDF time stepping, incremental Chorin projection, finite elements, and IEQ with mass lumping.The scheme is energy stable and preserves the symmetry and trace-free structure of Q.Source: §1 · pp. 3–4
• Convergence analysis must reconcile two projection velocities with different bounds and establish strong convergence of ∇Q without a discrete H^2 bound.The proof uses the structure of the equation for H to obtain precompactness of ∇Q; the velocity approximations are combined using tools from [27] [54] [53].Source: §1 · pp. 3–4
• IEQ's field-valued auxiliary variable enters the energy balance locally, which the convergence argument exploits; mass lumping also enables the discrete energy inequality and exact cancellation of terms.The local energy-balance role distinguishes the chosen IEQ approach from SAV's scalar auxiliary variable and nonlocal coupling.Source: §1 · p. 4
• The scheme reproduces smooth-solution spatial and temporal accuracy and qualitatively captures +1 defect splitting and skyrmion transport and deformation under a constant pressure gradient.The skyrmion simulation uses a constant pressure gradient to induce the flow.Source: §1 · p. 4

2. Preliminaries

• The preliminaries introduce notation and the weak-solution framework used for the Beris–Edwards system.The weak-solution notion is motivated by Leray–Hopf solutions for incompressible Navier–Stokes equations.Source: §2 · p. 4, §2.2 · p. 5

2.1. Notation.

• L^p(0, T; X) denotes strongly measurable X-valued functions that are L^p-integrable in time; the text also fixes norm, duality, and inner-product notation.The Leray projector maps L^2 vector fields to their divergence-free component.Source: §2.1 · p. 4

2.2. Weak solutions.

• A weak solution has symmetric, trace-free Q and H, divergence-free velocity, satisfies the distributional system and energy inequality, and is right-continuous at zero in the stated norms.This notion does not require Q ∈ L^2(0, T; H^2(Ω)), unlike the notions in.Source: §2.2 · p. 5

3.2. Spatial discretization.

• The spatial scheme uses conforming finite elements on a shape-regular, quasi-uniform mesh and mass lumping; its projection step is a Poisson problem requiring no inf-sup condition.The pressure is computed from the Poisson formulation, and the final velocity follows from the projection equation.Source: §3 · p. 5, §3.2 · pp. 6, 9

3.1. Time discretization.

• The time discretization uses formally second-order BDF2 with Chorin’s two-step projection, modifying s because the intermediate velocity need not be divergence free.When the velocity is divergence free, the modified s agrees with S.Source: §3.1 · p. 6

3.3. Solvability and properties of the scheme.

• Each time step is uniquely solvable with no relation between the time step and mesh size required for solvability or stability; convergence analysis requires h^2 = o(Δt).Source: §3.3 · p. 10
• The scheme propagates symmetry and trace-free structure: symmetric, trace-free initial Q yields symmetric, trace-free Q and molecular field H at subsequent steps.Source: §3.3 · pp. 11–12
• If the initial Q-tensor is symmetric and trace-free, the scheme preserves these properties for subsequent Q-tensors.The proof also establishes that any solution of the prediction step has a symmetric, trace-free Q-tensor.Source: §3.3 · pp. 13–14
• The scheme has a unique solution for each time step, with the prediction and projection steps well-posed.The projection step also yields a pressure satisfying the weak divergence constraint.Source: §3.3 · pp. 13–15

3.4. Energy stability.

• The approximations computed by the scheme satisfy a discrete energy estimate.The estimate is used to obtain a priori bounds and pre-compactness of the approximations.Source: §3.4 · p. 15

• The energy estimates yield weakly convergent subsequences, and the weak velocity limit is divergence free.The divergence-free conclusion uses that the limit lies in L2(0,T;H1_0(Ω)).Source: §4 · pp. 20–21
• The velocity approximations share the same limit, allowing time equicontinuity for one approximation to combine with spatial regularity for another.The Q-tensor approximations Qh and Q̄h also have the same limit; the argument addresses the fact that the available velocity estimates apply to different approximations.Source: §4 · pp. 20–21
• Under h^σ ≤ CΔt for 2 ≤ σ ≤ k+1, the velocity and Q-tensor approximations are weakly equicontinuous in time up to errors vanishing as h, Δt → 0.The same equicontinuity conclusion applies to rh when tested against H2(Ω) functions.Source: §4 · p. 26

Convergence analysis

• As h and Δt tend to zero with h^2 = o(Δt), the scheme’s approximations converge along a subsequence to a weak solution of the Beris–Edwards system.The converging sequences include the velocity approximations, Q-tensor approximations, molecular field, and auxiliary variable.Source: §4 · p. 29
• The limit velocity, Q-tensor, and auxiliary variable are weakly continuous in time with values in L2(Ω).The weak continuity of the auxiliary-variable limit is also stated for every t ∈ [0,T].Source: §4 · p. 29
• Convergence of the L2 norm of ∇Qh, combined with weak convergence, yields strong convergence in L2([0, T] × Ω) up to a subsequence; Lemma 4.1 also gives strong convergence of Q̄h in L2(0, T; H1(Ω)).Source: §4 · p. 34
• The limits satisfy the distributional Beris–Edwards equations, with the available time regularity ∂tQ ∈ L2(0,T; L6/5(Ω)) sufficient for the equivalence argument.The proof establishes the distributional formulation through the limiting equations and states that the spatial regularity follows from the discrete energy bounds.Source: §4 · p. 37

4. Convergence analysis of the scheme

• The limit satisfies the energy inequality (2.3).This follows after identifying r = r(Q).Source: §4 · p. 40
• Continuity at zero follows from weak time continuity of u, r, and ∇Q, together with weak lower semicontinuity and the energy inequality.For ∇Q, the argument uses ∂tQ ∈ L2([0, T]; L6/5(Ω)) and ∇Q ∈ L∞(0, T; L2(Ω)).Source: §4 · pp. 40–41
• The limit Q and H are symmetric and trace-free almost everywhere, so (u, Q, H) is a weak solution of (1.2).The discrete tensors’ symmetry and trace-free structure is passed to the limit using the self-adjointness of Π.Source: §4 · p. 41

5. Computational results

• The experiments use lowest-order Taylor–Hood velocity-pressure elements and lowest-order Lagrange spaces for tensor- and scalar-valued variables.The implementation uses Firedrake; discrete Q symmetry is strongly enforced.Source: §5 · p. 41

5.1. Spatial convergence test.

• Table 1 reports L2 spatial convergence rates at final time T = 1, with results suggesting quadratic convergence.The test uses manufactured solutions on a uniform triangulation with Δt = 1/1000.Source: §5.1 · pp. 41–42, Table 1 · p. 42

5.2. Temporal convergence test.

• The L2 temporal rates are quadratic for p, Q, and H, while u and r converge at rates ≈1.8 and ≈2.8, respectively.Rates are reported at final time T = 1.Source: Table 2 · p. 42

5.3. Simulation: defect dynamics.

• The defect simulation captures an initial +1 point defect splitting into two +1/2 defects, with higher velocity magnitude near the resulting defects.The findings qualitatively agree with [60, Section 4.2].Source: §5.3 · p. 42, Figure 1 · p. 43, Figure 2 · p. 44

5.4. 2.5D simulation: skyrmion dynamics.

• In the 2.5D simulation, the skyrmion is transported left to right and deformed by flow, in qualitative agreement with [21].The setup uses a constant pressure gradient to induce flow; the model’s analysis does not immediately apply to this hybrid setting.Source: §5.4 · p. 43

Appendix A. Discrete Grönwall inequality

• The lemma is identified as useful for proving an energy bound when the fluid equation has a nonzero external source.It is cited from [25, Proposition 4.1].Source: Appendix A · p. 45
• The lemma concerns nonnegative sequences satisfying a recurrence, with parameters ν and Δt subject to 1 − νΔt > 0.The recurrence itself is stated across the lemma’s displayed conditions.Source: Appendix A · p. 45
• An additional conclusion is stated for the case in which the sequence {b_n} is monotonically increasing.The cited lemma gives this as a further implication.Source: Appendix A · p. 45

Appendix B. Properties of the mass-lumped inner product

• For finite element functions, the mass-lumped inner product is bounded by products of discrete Lp norms, which in turn are controlled by the corresponding continuous Lp norms.The stated bound applies for conjugate exponents p and p′.Source: Appendix B · p. 46
• The difference between the mass-lumped and standard inner products is estimated by an O(h^(1+ℓ)) term involving the gradient norm of f.The supplied estimate is stated for finite element functions f and g.Source: Appendix B · pp. 45–46
• For Lipschitz-continuous F, the L1 interpolation error of v_h F(w_h) is bounded by C h times Lp and gradient-Lp′ norms.The estimate also extends to vector- and tensor-valued functions with the appropriate contraction.Source: Appendix B · p. 48

Appendix C. A technical consequence following from the Rellich theorem

• The lemma asserts that finitely many L2 functions can approximate any L2 function to an arbitrary prescribed tolerance in the stated estimate.The tolerance is denoted η > 0.Source: Appendix C · p. 48
• The construction uses an L2-orthonormal basis of Dirichlet Laplacian eigenfunctions whose eigenvalues tend to infinity.The eigenfunctions can be chosen in H1_0(Ω) ∩ C∞(Ω).Source: Appendix C · p. 48
• Choosing N so that (1 + λ_N)^−1 < η^2 controls the tail after truncating the eigenfunction expansion.The passage presents this choice as part of the finite-dimensional approximation argument.Source: Appendix C · p. 49