Modified Characteristic Finite Element Method with Second-Order Spatial Accuracy for Solving Convection-Dominated Problem on Surfaces

1. Introduction

The convection-reaction-diffusion (CRD) equation holds significant importance in the field of fluid mechanics as it serves as a model that captures the interconnected processes of convection, reaction, and diffusion. This equation proves instrumental in elucidating various natural phenomena such as alterations in liquid pollutants upon discharge into rivers, the durability of reinforced concrete in seawater [1], and heat conduction, among others. However, certain practical phenomena frequently manifest in irregular domains, including biological films [2,3], pattern formation [4,5], surfactant transportation [6,7], evolution of colonies on irregular surfaces [8], and cell migration [9,10]. The governing equation of these phenomena is the CRD equation on surfaces, thus necessitating the exploration of numerical methods for solving these equations on surfaces.

The numerical methods for solving partial differential equations (PDEs) on surfaces can be divided into two main categories: mesh-free methods [11,12,13,14] and mesh-based methods [7,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. For mesh-free methods, the implementation of this particular method is relatively straightforward. However, there are challenges in stability and error estimation. Conversely, mesh-based methods are highly dependent on the mesh, but improve the stability of the numerical method. We focus here on the finite element method which is one of mesh-based methods for solving PDEs on surfaces. The mesh generation include two prevalent strategies: embedding the surface in the narrow-band domain [7,18,19,20,21,22,23,24] and directly discretizing the surface [26,27,28,29,30,31,32,33].

When convection is dominant, the effectiveness of classical finite element methods is reduced, yielding non-physical oscillations in the numerical solution. To alleviate these oscillations, various finite element methods have been introduced, such as Olshanskii et al. [24] have presented error estimates for a trace finite element method with streamline-upwind/Petrov-Galerkin stabilization. This represents the first residual-based stabilization methods for convection-diffusion equations on surfaces. However, the accuracy and stability of this method depend on the careful selection of the stabilisation parameters, which requires rigorous theoretical considerations. For the same problem, Simon and Tobiska provide a priori error estimation for a fitted finite element local projection stabilization which is symmetric stability terms in their study [28]. Recently, Xiao et al. developed gradient recovery-based adaptive techniques in [29,30]. Jin et al. [31] subsequently utilized their work [29,30] to address the convection-dominated problem on surfaces using mixed finite element methods and provided theoretical analysis.

The characteristic finite element method [7,32,34] is easier to operate than the above methods [24,28,29,30,31]. To implement this method, it is necessary to convert the CRD equation into a reaction-diffusion equation and interpolate the solution in the characteristic direction [34]. Unlike the regular domains [34], the point in the characteristic direction pass through the surface only once during a short time. Consequently, the solution situated beyond the surface in the characteristic direction needs to be reconstructed based on available information. This task is straightforward for the characteristic finite element method based on volume mesh [7]. However, it presents challenges for the face mesh-based method in Section 4 of [26]. In response, Xiao et al. [32] introduced a characteristic finite element method (CFEM) that relies on the face mesh. By incorporating mass lumping, the CFEM ensures the preservation of positivity, albeit at the cost of diminished accuracy. Here, we will propose a modified characteristic finite element method with second-order spatial accuracy for solving the CRD equation on surfaces. The temporal discretization strategy employed is the backward-Euler method, while the spatial mesh adopts the face mesh. To reconstruct the solution beyond the surface in the characteristic direction, we consider the Taylor expansion. However, the high accuracy is concomitant with a reduction in the stability. This will result in the transformation of our characteristic finite element method into an explicit-implicit method, thereby compromising its ability to maintain positivity. Then, the reasons for the deteriorate in spatial accuracy of the characteristic finite element method based on face mesh are analyzed.

The rest of this paper is organized as follows. We start in Section 2 with an introduction to differential operators, Green’s theorem on surface and the CRD equation on surfaces. The reaction-diffusion equation with a characteristic directional derivative is then given. Additionally, we propose a modified characteristic finite element method based on Taylor expansion in Section 3. Then the reconstruction methods proposed by us and the CFEM are examined. Lastly, in Section 4, we present a series of numerical examples to assess the disparities between our scheme and the CFEM. The conclusion of finding is provided in Section 5.

2. Preliminaries

This section aims to present a thorough introduction to surface operators, Green’s theorem on surface and the CRD equation involved. Subsequently, the CRD equation is converted into a reaction-diffusion equation with a characteristic directional derivative employing a characteristic finite element method.

2.1. Surface Operators and Green’s Theorem on Surface

Let $\Gamma =\{\mathbf{x}\in {\mathbb{R}}^3|\psi \left(\mathbf{x}\right)=0\}$ be a connected and oriented surface without boundary for such $\psi \in C^2\left({\mathbb{R}}^3\right)$. For $\mathbf{x}\in \Gamma$, $P\left(\mathbf{x}\right)=I-\mathbf{n}\left(\mathbf{x}\right){\mathbf{n}}^{\mathrm{T}}\left(\mathbf{x}\right)$ denotes the tangential projection operator where $\mathbf{n}\left(\mathbf{x}\right)=\frac{\nabla \psi \left(\mathbf{x}\right)}{|\nabla \psi (\mathbf{x}\left)\right|}$ in the unit vector normal to $\Gamma$. The tangential gradient of $f\in C^2(\Gamma )$ is given by

$${\nabla }_{\Gamma }f\left(\mathbf{x}\right)=P\left(\mathbf{x}\right)·\nabla f\left(\mathbf{x}\right)\mathbf{x}\in \Gamma ,$$

(1)

and the Laplace-Beltremi operator is defined as

$${\Delta }_{\Gamma }f\left(\mathbf{x}\right)={\nabla }_{\Gamma }·{\nabla }_{\Gamma }f\left(\mathbf{x}\right)f\left(\mathbf{x}\right)\in C^2(\Gamma ).$$

(2)

Let H be the mean curvature of $\Gamma$, Green’s Theorem on any hypersurface is as follows.

Theorem 1. 

[26] If the boundary $\partial \Gamma$ of a hypersurface $\Gamma \subset {\mathbb{R}}^{n+1},n=1,2,...,$ is smooth, then $f\in C^1\left(\overline{\Gamma }\right)$ satisfy that

$${\int }_{\Gamma }{\nabla }_{\Gamma }fdA={\int }_{\Gamma }Hf\mathbf{n}dA+{\int }_{\partial \Gamma }f\mathbf{\nu }d\sigma ,$$

where $\mathbf{\nu }$ is the co-normal of Γ, $dA$ and $d\sigma$ are the measure of Γ and $\partial \Gamma$, respectively.

By choosing the inner product ${(u,v)}_{\Gamma }={\int }_{\Gamma }uv\mathrm{d}A$, the norms $\left|\right|·{\left|\right|}_{L^2(\Gamma )}$ and $\left|\right|·{\left|\right|}_{H^1(\Gamma )}$ are defined by

$${\left|\right|u\left|\right|}_{L^2(\Gamma )}=\sqrt{{(u,u)}_{\Gamma }}{\left|\right|u\left|\right|}_{H^1(\Gamma )}^2={\left|\right|u\left|\right|}_{L^2(\Gamma )}^2+\left|\right|{\nabla }_{\Gamma }u{\left|\right|}_{L^2(\Gamma )}^2.$$

(3)

2.2. The Convection-Reaction-Diffusion Equations on Surface

In this paper, we will focus on the following equation

$$\left\{\begin{array}{cc}{\partial }_{t}u+\mathit{\beta }·{\nabla }_{\Gamma }u-ϵ{\Delta }_{\Gamma }u+\mu u=f\hfill & (\mathbf{x},t)\in \Gamma \times (0,T],\hfill \\ u(\mathbf{x},0)=u_0\left(\mathbf{x}\right)\hfill & \mathbf{x}\in \Gamma ,\hfill \end{array}\right.$$

(4)

where the diffusion parameter $ϵ$ and the reaction coefficient $\mu$ are positive constants. The convection velocity $\mathit{\beta }$ assumed to be time-independent and continuous function satisfying

$$\mu -\frac{1}{2}{\nabla }_{\Gamma }·{\mathit{\beta }}_{\Gamma }\ge \kappa >0.$$

(5)

According to orthogonality, the convection term in (4) can be modified as

$$\begin{array}{cc}\hfill \mathit{\beta }·{\nabla }_{\Gamma }u& ={\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }u={\mathit{\beta }}_{\Gamma }·(P\nabla )u={\mathit{\beta }}_{\Gamma }·(I-P+P)\nabla u={\mathit{\beta }}_{\Gamma }·\nabla u,\hfill \end{array}$$

(6)

where ${\mathit{\beta }}_{\Gamma }$ represents the projection of $\mathit{\beta }$ onto the tangent plane of $\Gamma$. This means that we only need to pay attention to the tangential component ${\mathit{\beta }}_{\Gamma }$ of $\mathit{\beta }$ at the surface $\Gamma$.

2.3. The Reaction-Diffusion Equation with Characteristic Directional Derivative

For a small positive parameter $\delta$, let $U(\Gamma ,\delta )$ be the neighbourhood of the surface $\Gamma$ and $t^n=n\Delta t,n=0,1,...,N$, parameter $\Delta t=T/N$ is the time step with $\Delta t\le \delta$. We extend $u(\mathbf{x},t^n)$ from surface $\Gamma$ to neighborhood $U(\Gamma ,\delta )$. Assume there is a characteristic direction $\tau$ in $U(\Gamma ,\delta )\times [t^{n-1},t^n]$ and the point $\left(\mathbf{\chi }\right(t),t)$ on direction $\tau$ satisfy that

$$\left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\chi }\left(t\right)={\mathit{\beta }}_{\Gamma }t\in [t^{n-1},t^n),}\hfill \\ \mathbf{x}=\mathbf{\chi }\left(t^n\right).\hfill \end{array}\right.$$

(7)

Integrating the characteristic equation (7) over the interval from $t^{n-1}$ to $t^n$, the backtracking characteristic point

$$\mathbf{\chi }\left(t^{n-1}\right)=\mathbf{x}-\Delta t{\mathit{\beta }}_{\Gamma }\equiv \check{\mathbf{x}}.$$

(8)

Thanks to (7), the characteristic directional derivative at point $\left(\mathbf{\chi }\right(t),t)$ is

$$\begin{array}{c}\hfill {\partial }_{\tau }u=\frac{\mathrm{d}}{\mathrm{d}t}u(\mathbf{\chi }\left(t\right),t)={\partial }_{t}u+\nabla u·\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\chi }\left(t\right)={\partial }_{t}u+{\mathit{\beta }}_{\Gamma }·\nabla u.\end{array}$$

(9)

It follows from (4) and (9) that

$$\left\{\begin{array}{cc}{\partial }_{\tau }u-\epsilon {\Delta }_{\Gamma }u+\mu u=f\hfill & (\mathbf{x},t)\in \Gamma \times (0,T],\hfill \\ u(\mathbf{x},0)=u_0\left(\mathbf{x}\right)\hfill & \mathbf{x}\in \Gamma .\hfill \end{array}\right.$$

(10)

Several discrete schemes of the characteristic directional derivative ${\partial }_{\tau }u\left(t^n\right)$ are known, including the backward-Euler method, the Crank-Nicolson method, and the Runge-Kutta method. In this paper we employ the backward-Euler scheme:

$${\partial }_{\tau }u(\mathbf{x},t^n)=\frac{u(\mathbf{x},t^n)-u(\check{\mathbf{x}},t^{n-1})}{\Delta t}+O(\Delta t).$$

(11)

The backtracking characteristic point $\check{\mathbf{x}}$ is evidently not situated on $\Gamma$. Consequently, it is of utmost significance to reconstruct $u(\check{\mathbf{x}},t^{n-1})$ by utilizing the available information on surfaces.

3. A Modified Characteristic Finite Element Method (MCFEM) Based on Taylor Expansion

In this section, we will introduce a modified characteristic finite element method (MCFEM) which employs the Taylor expansion to obtain the solution beyond the surface in the characteristic direction. Then, we will provide a specific discrete formulation of the MCFEM and analysis its stability.

3.1. The Reconstruction Method Based on Taylor Expansion

The approximation of the solution in the characteristic direction is heavily mesh-dependent. If the volume mesh [7,18,19,20,21,22,23,24] is used, the solution beyond the surface can be reconstructed using interpolation. Hovever, the face mesh in Section 4 of [26] cannot generate a narrow band which contains surfaces similar to the volume mesh. Consequently, if a point in the characteristic direction is situated on the current mesh, it must have resided beyond the mesh at the previous moment. This presents novel challenges to the reconstruction method.

At present, the reconstruction method based on face mesh is only involved by Xiao et al. [32], as illustrated in Figure 1. They will identify the nearest mesh point ${\mathbf{x}}_c$ and project $\check{\mathbf{x}}$ vertically onto the tangent plane ${\mathrm{T}}_c$ of ${\mathbf{x}}_c$ as ${\mathbf{x}}^{*}$. Subsequently, the mesh points on each element containing ${\mathbf{x}}_c$ are extended into the tangent plane ${\mathrm{T}}_c$ along their normal direction to form a new local element. Linear interpolation is then employed to approximate $u\left({\mathbf{x}}^{*}\right)$ within these local elements, which is used to approximate $u\left(\check{\mathbf{x}}\right)$. It is evident that the CFEM in [32] fails to account for the discrepancy between $u\left({\mathbf{x}}^{*}\right)$ and $u\left(\check{\mathbf{x}}\right)$, resulting in a spatial accuracy lower than the second-order.

We will suggest a more accurate method for reconstruction the solution beyond the surface. It should be noted that the backtracking characteristic point $\check{\mathbf{x}}$ can be situated within the tangent plane ${\mathrm{T}}_{\mathbf{x}}$ of point $\mathbf{x}\in \Gamma$ through the selection of the suitable characteristic direction $\tau$. Considering the Taylor expansion within the tangent plane ${\mathrm{T}}_{\mathbf{x}}$ and (8), we have

$$\begin{array}{cc}\hfill u(\check{\mathbf{x}},t)& =u(\mathbf{x},t)+(\check{\mathbf{x}}-\mathbf{x})·\nabla u(\mathbf{x},t)+O\left(\right|\check{\mathbf{x}}-\mathbf{x}{|}^2)\hfill \\ \hfill & =u(\mathbf{x},t)-\Delta t{\mathit{\beta }}_{\Gamma }\left(\mathbf{x}\right)·\nabla u(\mathbf{x},t)+O(\Delta t^2)\hfill \\ \hfill & =\check{\tilde{u}}(\mathbf{x},t)+O(\Delta t^2),\hfill \end{array}$$

(12)

where $\check{\tilde{u}}\left(\mathbf{x}\right)$ is an approximation of $u\left(\check{\mathbf{x}}\right)$ with second-order accuracy in time.

Taking this approximation, we obtain the reconstruction method based on Taylor expansion, as shown in Figure 2. First, the backtracking characteristic point $(\check{\mathbf{x}},t^{n-1})$ of point $(\mathbf{x},t^n)$ on $\Gamma$ is found along the characteristic direction $\tau$. Due to (8), it is easy to know that the backtracking characteristic point is given by $\check{\mathbf{x}}=\mathbf{x}-\Delta t{\mathit{\beta }}_{\Gamma }\left(\mathbf{x}\right)$. Second, the reconstructed function $\check{\tilde{u}}(\mathbf{x},t^{n-1})=u(\mathbf{x},t^{n-1})-\Delta t{\mathit{\beta }}_{\Gamma }\left(\mathbf{x}\right)·\nabla u(\mathbf{x},t^{n-1})$ is obtained by using the Taylor expansion of $u(\check{\mathbf{x}},t^{n-1})$ at point $(\mathbf{x},t^{n-1})$.

3.2. Temporal Discretization of the MCFEM

Let ${\overline{\partial }}_{\tau }u\left(t^n\right)=\frac{u\left(t^n\right)-\check{\tilde{u}}\left(t^{n-1}\right)}{\Delta t}$ be an approximation of the characteristic directional derivative ${\partial }_{\tau }u\left(t^n\right)$. An estimate between ${\partial }_{\tau }u\left(t^n\right)$ and ${\overline{\partial }}_{\tau }u\left(t^n\right)$ is provided below.

Theorem 2. 

If $u\in C^2(U(\Gamma ,\delta )\times [0,T])$ satisfies equation (7), then the following estimate holds:

$$\begin{array}{c}\hfill \left|\right|{\partial }_{\tau }u\left(t^n\right)-{\overline{\partial }}_{\tau }u\left(t^n\right){\left|\right|}_{L^2(\Gamma )}\le {\int }_{t^{n-1}}^{t^n}\left|\right|{\partial }_{tt}{u\left|\right|}_{L^2(\Gamma )}+\left|\right|{\partial }_t\left({\partial }_{\tau }u\right){\left|\right|}_{L^2(\Gamma )}\mathrm{d}t,1\le n\le N.\end{array}$$

(13)

Proof. 

By the definition of operator ${\partial }_{\tau }u\left(t^n\right)$ and ${\overline{\partial }}_{\tau }u\left(t^n\right)$

$$\begin{array}{cc}\hfill \left|\right|{\partial }_{\tau }u\left(t^n\right)-{\overline{\partial }}_{\tau }u\left(t^n\right){\left|\right|}_{L^2(\Gamma )}& =\frac{1}{\Delta t}\left|\right|\Delta t({\partial }_{t}u\left(t^n\right)+{\mathit{\beta }}_{\Gamma }·\nabla u\left(t^n\right))-(u\left(t^n\right)-u\left(t^{n-1}\right)+\Delta t{\mathit{\beta }}_{\Gamma }·\nabla u\left(t^{n-1}\right)){\left|\right|}_{L^2(\Gamma )}\hfill \\ & =\frac{1}{\Delta t}\left|\right|\left[\Delta t{\partial }_{t}u\left(t^n\right)-(u\left(t^n\right)-u\left(t^{n-1}\right))\right]+\left[\Delta t{\mathit{\beta }}_{\Gamma }·\nabla (u\left(t^n\right)-u\left(t^{n-1}\right))\right]{\left|\right|}_{L^2(\Gamma )}\hfill \\ & :=\frac{1}{\Delta t}\left|\right|{\mathrm{T}}_1+{\mathrm{T}}_2{\left|\right|}_{L^2(\Gamma )}.\hfill \end{array}$$

(14)

For ${\mathrm{T}}_1$, we see that

$$\begin{array}{cc}\hfill {\mathrm{T}}_1& =\Delta t{\partial }_{t}u\left(t^n\right)-(u\left(t^n\right)-u\left(t^{n-1}\right))\hfill \\ \hfill & ={\int }_{t^{n-1}}^{t^n}{\partial }_{t}u\left(t^n\right)-{\partial }_{t}u\left(t\right)dt\hfill \\ \hfill & =(t-t^{n-1})({\partial }_{t}u\left(t^n\right)-{\partial }_{t}u\left(t\right)){|}_{t^{n-1}}^{t^n}+{\int }_{t^{n-1}}^{t^n}(t-t^{n-1}){\partial }_{tt}u\left(t\right)dt\hfill \\ \hfill & ={\int }_{t^{n-1}}^{t^n}(t-t^{n-1}){\partial }_{tt}u\mathrm{d}t.\hfill \end{array}$$

(15)

For ${\mathrm{T}}_2$, we have

$$\begin{array}{cc}\hfill {\mathrm{T}}_2& =\Delta t{\mathit{\beta }}_{\Gamma }·\nabla (u\left(t^n\right)-u\left(t^{n-1}\right))\hfill \\ \hfill & =\Delta t{\mathit{\beta }}_{\Gamma }·\nabla \left({\int }_{t^{n-1}}^{t^n}{\partial }_{t}u\mathrm{d}t\right)\hfill \\ \hfill & =\Delta t{\int }_{t^{n-1}}^{t^n}{\mathit{\beta }}_{\Gamma }·\nabla \left({\partial }_{t}u\right)\mathrm{d}t.\hfill \end{array}$$

(16)

By substituting (15) and (16) into (14), the bound of $|{\partial }_{\tau }u\left(t^n\right)-{\overline{\partial }}_{\tau }u\left(t^n\right)|$ can be obtained by

$$\begin{array}{cc}\hfill \left|\right|{\partial }_{\tau }u\left(t^n\right)-\frac{u\left(t^n\right)-\check{\tilde{u}}\left(t^{n-1}\right)}{\Delta t}{\left|\right|}_{L^2(\Gamma )}& =\frac{1}{\Delta t}\left|\right|{\int }_{t^{n-1}}^{t^n}(t-t^{n-1}){\partial }_{tt}u+\Delta t{\mathit{\beta }}_{\Gamma }·\nabla \left({\partial }_{t}u\right)\mathrm{d}t{\left|\right|}_{L^2(\Gamma )}\hfill \\ \hfill & =\frac{1}{\Delta t}\left|\right|{\int }_{t^{n-1}}^{t^n}(t-t^n){\partial }_{tt}u+\Delta t({\partial }_{tt}u+{\mathit{\beta }}_{\Gamma }·\nabla \left({\partial }_{t}u\right))\mathrm{d}t\left|\right|\hfill \\ \hfill & =\frac{1}{\Delta t}\left|\right|{\int }_{t^{n-1}}^{t^n}(t-t^n){\partial }_{tt}u+\Delta t{\partial }_t\left({\partial }_{\tau }u\right)\mathrm{d}t{\left|\right|}_{L^2(\Gamma )}\hfill \\ \hfill & \le {\int }_{t^{n-1}}^{t^n}\left|\right|{\partial }_{tt}u\left|\right|+\left|\right|{\partial }_t\left({\partial }_{\tau }u\right){\left|\right|}_{L^2(\Gamma )}\mathrm{d}t.\hfill \end{array}$$

(17)

This complete the proof. □

Let $u^n$ is an approximation of $u\left(t^n\right)$ and ${\check{\tilde{u}}}^{n-1}=u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·\nabla u^{n-1}$. The temporal discretization of (10) is to find $u^n\in H^1(\Gamma ),n=1,2,...,N$, recursively

$$\left\{\begin{array}{c}{(\overline{{\partial }_{\tau }}{u}^n,v)}_{\Gamma }+ϵ{({\nabla }_{\Gamma }{u}^n,{\nabla }_{\Gamma }v)}_{\Gamma }+{(\mu u^n,v)}_{\Gamma }={(f^n,v)}_{\Gamma }\forall v\in H^1(\Gamma ),\hfill \\ u^0=u_0,\hfill \end{array}\right.$$

(18)

where $f^n=\frac{1}{\Delta t}{\int }_{t^{n-1}}^{t^n}f\left(t\right)\mathrm{d}t$. The stability of the problem (18) is demonstrated in the following manner.

Theorem 3. 

Suppose that ${\left\{u^n\right\}}_{n=1,2,...,N}\subset H^1(\Gamma )$ satisfy (18). If $\Delta t\le \frac{2ϵ}{\left|\right|{\mathit{\beta }}_{\Gamma }{\left|\right|}_{L^2(\Gamma )}^2}$, then the following inequality holds:

$$\begin{array}{cc}\hfill max\left\{\left|\right|u^n-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2,\left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2+2\Delta t\mu \left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2+2\Delta tϵ\left|\right|{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\right\}& \le e^T\left(\left|\right|u^0-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^0{\left|\right|}_{L^2(\Gamma )}^2\right.\hfill \\ & \left.+{\int }_0^T{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^{2}dt\right),\hfill \end{array}$$

(19)

for all n.

Proof. 

By the definition of $\overline{{\partial }_{\tau }}{u}^n$, (18) can be written as

$$\begin{array}{cc}\hfill {(\frac{u^n-[u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}]}{\Delta t},v)}_{\Gamma }+ϵ{({\nabla }_{\Gamma }{u}^n,{\nabla }_{\Gamma }v)}_{\Gamma }+\mu {(u^n,v)}_{\Gamma }& ={(f^n,v)}_{\Gamma }.\hfill \end{array}$$

(20)

Taking $v=2\Delta t{u}^n$ into (20), we have

$$\begin{array}{cc}\hfill {(2u^n+2\Delta t\mu u^n,u^n)}_{\Gamma }+2\Delta tϵ\left|\right|{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2& =2{(u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}+\Delta t{f}^n,u^n)}_{\Gamma }\hfill \\ & \le \left|\right|u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}+\Delta t{f}^n{\left|\right|}_{L^2(\Gamma )}^2+\left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2\hfill \\ & \le \left|\right|u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}{\left|\right|}_{L^2(\Gamma )}^2+2\Delta t{(u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1},f^n)}_{\Gamma }\hfill \\ & +\Delta t^2\left|\right|f^n{\left|\right|}_{L^2(\Gamma )}^2+\left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2\hfill \\ & \le (1+\Delta t)\left|\right|u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}{\left|\right|}_{L^2(\Gamma )}^2+\Delta t(1+\Delta t)\left|\right|f^n{\left|\right|}_{L^2(\Gamma )}^2+\left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2.\hfill \end{array}$$

(21)

For $\left|\right|f^n{\left|\right|}_{L^2(\Gamma )}^2$, we see that

$$\begin{array}{cc}\hfill \left|\right|f^n{\left|\right|}_{L^2(\Gamma )}^2& ={\int }_{\Gamma }{\left(\frac{1}{\Delta t}{\int }_{t^{n-1}}^{t^n}f\left(t\right)dt\right)}^2\mathrm{d}A\hfill \\ \hfill & =\frac{1}{\Delta t^2}{\int }_{\Gamma }{\int }_{t^{n-1}}^{t^n}{\int }_{t^{n-1}}^{t^n}f\left(t\right)f\left(s\right)\mathrm{d}s\mathrm{d}t\mathrm{d}A\hfill \\ \hfill & \le \frac{1}{\Delta t^2}{\int }_{\Gamma }{\int }_{t^{n-1}}^{t^n}{\int }_{t^{n-1}}^{t^n}\frac{{\left(f\left(t\right)\right)}^2+{\left(f\left(s\right)\right)}^2}{2}\mathrm{d}s\mathrm{d}t\mathrm{d}A\hfill \\ \hfill & \le \frac{1}{\Delta t}{\int }_{t^{n-1}}^{t^n}{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^2\mathrm{d}t.\hfill \end{array}$$

(22)

Taking above inequality into (21), we get

$$\begin{array}{cc}\hfill \left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2+2\Delta t{(\mu u^n,u^n)}_{\Gamma }+2\Delta tϵ\left|\right|{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2& \le (1+\Delta t)\left(\left|\right|u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}{\left|\right|}_{L^2(\Gamma )}^2+{\int }_{t^{n-1}}^{t^n}{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^2\mathrm{d}t\right).\hfill \end{array}$$

(23)

It follows (5) that the bound

$$\begin{array}{c}\hfill -2{({\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n,u^n)}_{\Gamma }={({\nabla }_{\Gamma }·{\mathit{\beta }}_{\Gamma },{\left(u^n\right)}^2)}_{\Gamma }\le 2{(\mu u^n,u^n)}_{\Gamma },\end{array}$$

(24)

is hold. If $\Delta t\le \frac{2ϵ}{\left|\right|{\mathit{\beta }}_{\Gamma }{\left|\right|}_{L^2(\Gamma )}^2}$ holds, we have

$$\begin{array}{c}\hfill \Delta t\left|\right|{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\le \Delta t\left|\right|{\mathit{\beta }}_{\Gamma }{\left|\right|}_{L^2(\Gamma )}^2·\left|\right|{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\le 2ϵ\left|\right|{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2.\end{array}$$

(25)

Combining inequalities (24), (25) and (23), the key idea of proving Theorem 3 is given by

$$\begin{array}{cc}\hfill \left|\right|u^n-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2& ={(u^n,u^n)}_{\Gamma }-2{(\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n,u^n)}_{\Gamma }+\Delta t^2\left|\right|{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\hfill \\ \hfill & \le {(u^n,u^n)}_{\Gamma }+2\Delta t{(\mu u^n,u^n)}_{\Gamma }+2\Delta tϵ\left|\right|{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\hfill \\ \hfill & \le (1+\Delta t)\left(\left|\right|u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}{\left|\right|}_{L^2(\Gamma )}^2+{\int }_{t^{n-1}}^{t^n}{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^2\mathrm{d}t\right).\hfill \end{array}$$

(26)

After repeated application the inequality (26), we have

$$\begin{array}{cc}\hfill \left|\right|u^n-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2& \le {(1+\Delta t)}^n\left|\right|u^0-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^0{\left|\right|}_{L^2(\Gamma )}^2+\sum _{i=1}^n{\int }_{t^{i-1}}^{t^i}{(1+\Delta t)}^{n+1-i}{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^2\mathrm{d}t\hfill \\ & \le {(1+\Delta t)}^n\left(\left|\right|u^0-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^0{\left|\right|}_{L^2(\Gamma )}^2+{\int }_0^T{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^2\mathrm{d}t\right).\hfill \end{array}$$

(27)

Due to ${(1+\Delta t)}^n\le {(1+\Delta t)}^N={(1+\Delta t)}^{\frac{T}{\Delta t}}\le e^T$, the following inequality obviously

$$\begin{array}{c}\hfill \left|\right|u^n-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\le e^T\left(\left|\right|u^0-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^0{\left|\right|}_{L^2(\Gamma )}^2+{\int }_0^T{\left|\right|f\left(t\right)\left|\right|}_{L^2(\Gamma )}^2\mathrm{d}t\right),\end{array}$$

(28)

holds. By combining (28), (27) and (23), the proof of Theorem 3 is completed. □

3.3. The Surface Finite Element Method

Let ${\left\{{\Gamma }_h\right\}}_{h>0}$ be a family of discrete surfaces which composed of plane open triangles $K_j$ with edge $\partial K_j$ and vertexes ${\mathbf{x}}_l,l=1,...,N_h^v$. The point ${\mathbf{x}}_l$ is also the vertex points of curved open triangle $\mathcal{K}\in \Gamma$ such that

$$\underset{j=1}{\overset{N_h^{\mathcal{T}}}{\cup }}{\overline{\mathcal{K}}}_j=\Gamma ,$$

and for $j\ne k$, ${\overline{\mathcal{K}}}_j\cap {\overline{\mathcal{K}}}_k=\varnothing$ or common curved edge of ${\overline{\mathcal{K}}}_j$ and ${\overline{\mathcal{K}}}_k$ or common vertex of ${\overline{\mathcal{K}}}_j$ and ${\overline{\mathcal{K}}}_k$. For an interior edge $E_j,j=1,2,...,N_h^E$, there are two triangles $K_l^j$ and $K_r^j$ such that $\partial K_l^j\cap \partial K_r^j=E_j$. Following [26], we adopt the projection ${\mathcal{P}}_{{\Gamma }_h}:\Gamma \to {\Gamma }_h$ which is Lipschitz continuous and ${\mathcal{P}}_{{\Gamma }_h}\left({\mathcal{K}}_j\right)=K_j$ for any triangle element $K_j\subset {\Gamma }_h$. For any $f\in C^0(\Gamma )$, its projection on $\Gamma$ is given by $f_{{\mathcal{P}}_{{\Gamma }_h}}=f\circ {\mathcal{P}}_{{\Gamma }_h}^{-1}$. The projection of ${\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}$ on the tangent plane of the discrete surface ${\Gamma }_h$ is denoted as ${\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h},{\Gamma }_h}$. Let finite dimensional space $S_h$ be a continuous function space on ${\Gamma }_h$ that is linear on each triangle $K_j$. Considering the variational problem (18), we get the MCFEM of equation (10): For $n=1,2,...,N$, find $u_h^n=u_h(\mathbf{x},t^n)\in S_h$ such that

$$\left\{\begin{array}{cc}{(\overline{{\partial }_{\tau }}{u}_h^n,v_h)}_{{\Gamma }_h}+\epsilon {({\nabla }_{{\Gamma }_h}{u}_h^n,{\nabla }_{{\Gamma }_h}{v}_h)}_{{\Gamma }_h}+\mu {(u_h^n,v_h)}_{{\Gamma }_h}={(f_{{\mathcal{P}}_{{\Gamma }_h}}^n,v_h)}_{{\Gamma }_h}\hfill & \forall v_h\in S_h,\hfill \\ u_h^0\left(\mathbf{x}\right)=I_h{u}_0\left(\mathbf{x}\right),\hfill \end{array}\right.$$

(29)

where $\overline{{\partial }_{\tau }}{u}_h^n=\frac{u_h^n-[u_h^{n-1}-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h},{\Gamma }_h}·\nabla u_h^{n-1}]}{\Delta t}$, $f_{{\mathcal{P}}_{{\Gamma }_h}}^n=\frac{1}{\Delta t}{\int }_{t^{n-1}}^{t^n}{f}_{{\mathcal{P}}_{{\Gamma }_h}}\left(t\right)\mathrm{d}t$ and $I_h$ is a piecewise linear interpolation operator.

Although the method (29) is based on the idea of characteristic finite element, it utilizes Taylor’s expansion to restructure $u\left(\check{\mathbf{x}}\right)$. Consequently, the MCFEM (29) degenerates into an explicit-implicit method, imposing stringent restriction on the mesh size h and time step $\Delta t$. Before analyzing the stability of the MCFEM (29), we need to restrict the mesh size.

Theorem 4. 

For any point ${\mathbf{x}}_0\in \Gamma$, there are two triangles $K_l,K_r\subset {\Gamma }_h$ such that ${\mathbf{x}}_0\in \partial {\mathcal{K}}_l\cap \partial {\mathcal{K}}_r$ or ${\mathbf{x}}_0$ is not the vertex points of ${\mathcal{K}}_l$ and ${\mathcal{K}}_r$. ${\mathbf{\nu }}_r$ and ${\mathbf{\nu }}_l$ denote the unit outward normal vectors to $K_l$ and $K_r$, respectively. If the convection velocity ${\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}\in C^0\left({\Gamma }_h^3\right)$, then there exists a specific mesh size $h_{\kappa }$ such that the following inequality holds

$$\left|\left[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]\right]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)\right|\le \kappa \forall h<h_{\kappa },$$

where $\left[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h},{\Gamma }_h}]\right]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)=\underset{h\to 0^{+}}{lim}\left[{\mathbf{\nu }}_l·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_l}+{\mathbf{\nu }}_r·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_r}\right]$.

Proof. 

By the definition, we have

$$\begin{array}{cc}\hfill \left[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]\right]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)& =\underset{h\to 0^{+}}{lim}\left[{\mathbf{\nu }}_l·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_l}+{\mathbf{\nu }}_r·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_r}\right]\hfill \\ \hfill & =\underset{h\to 0^{+}}{lim}\left[{\mathbf{\nu }}_l·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_l}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_l}+{\mathbf{\nu }}_r·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_r}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_r}\right]\hfill \\ \hfill & =\underset{h\to 0^{+}}{lim}\left[({\mathbf{\nu }}_l·({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}-({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}·{\mathbf{n}}_{K_l}){\mathbf{n}}_{K_l})){|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_l}\right]\hfill \\ & +\underset{h\to 0^{+}}{lim}\left[({\mathbf{\nu }}_r·({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}-({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}·{\mathbf{n}}_{K_r}){\mathbf{n}}_{K_r})){|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_r}\right],\hfill \end{array}$$

(30)

where ${\mathbf{n}}_{K_l}$ and ${\mathbf{n}}_{K_r}$ are the unit vectors normal to $K_l$ and $K_r$, respectively. Considering the orthogonality, (30) can be rewritten as

$$\begin{array}{cc}\hfill \left[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]\right]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)& =\underset{h\to 0^{+}}{lim}\left[{\mathbf{\nu }}_l·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_l}+{\mathbf{\nu }}_r·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_r}\right].\hfill \end{array}$$

(31)

It is obvious that

$$\begin{array}{cc}\hfill \left[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]\right]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)& =\underset{h\to 0^{+}}{lim}\left[{\mathbf{\nu }}_l·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_l}+{\mathbf{\nu }}_r·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h}}{|}_{\mathbf{x}={\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)+h{\mathbf{\nu }}_r}\right]\hfill \\ \hfill & =\mathbf{\nu }·\mathit{\beta }{|}_{\mathbf{x}={\mathbf{x}}_0}-\mathbf{\nu }·\mathit{\beta }{|}_{\mathbf{x}={\mathbf{x}}_0}=0.\hfill \end{array}$$

(32)

Hence, the existence of the mesh size $h_{\kappa }$ required by Theorem 4 can be established. The proof is completed. □

Next, the stability of scheme (29) will be demonstrated.

Theorem 5. 

Assume that $\Delta t\le \frac{2ϵ}{\left|\right|{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2}$, $h\le h_{\kappa }$ and (5) hold. Then the solution $u_h^n$ of Problem (29) satisfies

$$\begin{array}{cc}\hfill max\left\{\left|\right|u_h^n-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2,\left|\right|u_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2+2\Delta t\mu \left|\right|u_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2+2\Delta tϵ\left|\right|{\nabla }_{{\Gamma }_h}{u}_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2\right\}& \le e^T\left(\left|\right|u_h^0-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^0{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2\right.\hfill \\ & \left.+{\int }_0^T\left|\right|f_{{\mathcal{P}}_{{\Gamma }_h}}\left(t\right){\left|\right|}_{L^2\left({\Gamma }_h\right)}^{2}dt\right),\hfill \end{array}$$

(33)

for $n=1,2,...,N.$

Proof. 

Choose $v_h=2\Delta t{u}_h^n$ in (29) so that

$$\begin{array}{cc}\hfill {(2u_h^n+2\Delta t·\mu u_h^n,u_h^n)}_{{\Gamma }_h}+2\Delta tϵ\left|\right|{\nabla }_{{\Gamma }_h}{u}_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2& =2{(u_h^{n-1}-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^{n-1}+\Delta t{f}_{{\mathcal{P}}_{{\Gamma }_h}}^n,u_h^n)}_{{\Gamma }_h}\hfill \\ & \le \left|\right|u_h^{n-1}-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^{n-1}+\Delta t{f}_{{\mathcal{P}}_{{\Gamma }_h}}^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2+\left|\right|u_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2.\hfill \end{array}$$

(34)

Similarly to Theorem 3, we have

$$\begin{array}{cc}\hfill \left|\right|u_h^n-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2+2\Delta t{(\mu u_h^n+{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}& \le (1+\Delta t)\left(\left|\right|u_h^{n-1}-\Delta t{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^{n-1}{\left|\right|}_{L^2\left({\Gamma }_h\right)}^2\right.\hfill \\ & +{\int }_{t^{n-1}}^{t^n}\left.\left|\right|f_{{\mathcal{P}}_{{\Gamma }_h}}\left(t\right){\left|\right|}_{L^2\left({\Gamma }_h\right)}^2\mathrm{d}t\right).\hfill \end{array}$$

(35)

We claim that ${(\mu u_h^n+{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}\ge 0$. Owing to (5), the following inequality holds

$$\begin{array}{cc}\hfill 2{(\mu u_h^n+{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}& \ge 2\kappa {(u_h^n,u_h^n)}_{{\Gamma }_h}+{({\nabla }_{{\Gamma }_h}·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}},{\left(u_h^n\right)}^2)}_{{\Gamma }_h}+2{({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}.\hfill \end{array}$$

(36)

Since the discrete surface ${\Gamma }_h$ is a piecewise linear approximation of $\Gamma$, the unit outward normal vectors to any adjacent triangles $K_l^j$ and $K_r^j$ are discontinuous at a common edge $E_j$. Consequently, the last two terms on the right side of inequality (36) are not equal to 0. Using Theorem 1, we obtain

$$\begin{array}{cc}\hfill {({\nabla }_{{\Gamma }_h}·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}},{\left(u_h^n\right)}^2)}_{{\Gamma }_h}+2{({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}& ={\int }_{{\Gamma }_h}{\nabla }_{{\Gamma }_h}\left({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}{\left(u_h^n\right)}^2\right)\mathrm{d}{A}_h\hfill \\ & =\sum _{j=1}^{N_h^{\mathcal{T}}}{\int }_{K_j}{\nabla }_{K_j}\left({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_j}}{\left(u_h^n\right)}^2\right)\mathrm{d}{A}_h\hfill \\ & =\sum _{j=1}^{N_h^{\mathcal{T}}}\left[{\int }_{K_j}{H}_{K_j}{\left(u_h^n\right)}^2{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_j}}·{\mathbf{n}}_{K_j}\mathrm{d}{A}_h+{\int }_{\partial K_j}{\left(u_h^n\right)}^2{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_j}}·{\mathbf{\nu }}_{\partial K_j}\mathrm{d}{\sigma }_h\right],\hfill \end{array}$$

(37)

where ${\mathbf{n}}_{K_j}$ is the unit vector normal to $K_j$, ${\mathbf{\nu }}_{\partial K_j}$ is unit normal vector of boundary $\partial K_j$ and tangent to $K_j$. Note that $\mathrm{d}{A}_h$ and $\mathrm{d}{\sigma }_h$ are the measure of ${\Gamma }_h$ and interior edge $\partial K_j$, respectively. The mean curvature of triangular element $K_j$ is defined as $H_{K_j}$. It follows form $K_j$ is a triangular plane that $H_{K_j}$ is equal to 0. By choosing the mesh size $h<h_{\kappa }$ and applying Theorem 4, a lower bound for the inequality (37) can be obtained by

$$\begin{array}{cc}\hfill {({\nabla }_{{\Gamma }_h}·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}},{\left(u_h^n\right)}^2)}_{{\Gamma }_h}+2{({\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}& =\sum _{j=1}^{N_h^{\mathcal{T}}}{\int }_{\partial K_j}{\left(u_h^n\right)}^2{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_j}}·{\mathbf{\nu }}_{\partial K_j}\mathrm{d}{\sigma }_h\hfill \\ & =\sum _{j=1}^{N_h^E}{\int }_{\partial K_l^j\cap \partial K_r^j}{\left(u_h^n\right)}^2\left[{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_l^j}}·{\mathbf{\nu }}_{\partial K_l^j}+{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,K_r^j}}·{\mathbf{\nu }}_{\partial K_r^j}\right]\mathrm{d}{\sigma }_h\hfill \\ & \ge -\kappa \sum _{j=1}^{N_h^E}{\int }_{\partial K_l^j\cap \partial K_r^j}{\left(u_h^n\right)}^2\mathrm{d}{\sigma }_h=-\kappa \sum _{j=1}^{N_h^{\mathcal{T}}}{\int }_{\partial K_j}{\left(u_h^n\right)}^2\mathrm{d}{\sigma }_h.\hfill \end{array}$$

(38)

If we plug inequality (38) into (36), we get

$$\begin{array}{cc}\hfill 2{(\mu u_h^n+{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}·{\nabla }_{{\Gamma }_h}{u}_h^n,u_h^n)}_{{\Gamma }_h}& \ge 2\kappa \sum _{j=1}^{N_h^{\mathcal{T}}}{\int }_{\partial K_j}{\left(u_h^n\right)}^2\mathrm{d}{\sigma }_h-\kappa \sum _{j=1}^{N_h^{\mathcal{T}}}{\int }_{\partial K_j}{\left(u_h^n\right)}^2\mathrm{d}{\sigma }_h\ge 0.\hfill \end{array}$$

(39)

Others are similar to Theorem 3. The proof of Theorem 5 is completed. □

3.4. The Analysis of Reconstruction Methods in MCFEM and CFEM

The reconstruction method proposed in Subsection 3.1 is different from that employed in the CFEM [32], as the scheme of the MCFEM does not involve projection. This method pre-processes the convection velocity $\mathit{\beta }$ to guarantee that the backtracking characteristic point $\check{\mathbf{x}}$ is situated on the tangent plane ${\mathrm{T}}_{\mathbf{x}}$ of the surface $\Gamma$. Taylor expansion is utilized in the MCFEM to maintain a second-order spatial accuracy. However, it also transforms the MCFEM into an explicit-implicit method, reducing its stability compared to the CFEM.

As described in Subsection 3.1, the reconstruction method employed by the CFEM involves projecting the backtracking characteristic point $\check{\mathbf{x}}$ and its closer mesh points onto tangent planes ${\mathrm{T}}_c$, thereby generating a novel local mesh. This also makes the CFEM more stable than the MCFEM. Since the characteristic finite element employs the definition domain of the solution u to be a neighbourhood $U(\Gamma ,\delta )$ containing the surface $\Gamma$, the value of u should be different at the same normal vector on $\Gamma$. This is contrary to the more widely accepted understanding in the use of conventional employs surface finite element method described in [26].

Employing the notations given in Subsection 3.1, the errors in the reconstruction method of the CFEM consist of $|u\left(\check{\mathbf{x}}\right)-u\left({\mathbf{x}}^{*}\right)|$ and $|u\left({\mathbf{x}}^{*}\right)-u_h\left({\mathbf{x}}^{*}\right)|$. Obviously,

$$|u\left(\check{\mathbf{x}}\right)-u\left({\mathbf{x}}^{*}\right)|\le \underset{\mathbf{x}\in U(\Gamma ,\delta )}{max}|\nabla u\left(\mathbf{x}\right)|·|\check{\mathbf{x}}-{\mathbf{x}}^{*}|.$$

(40)

Here, $|\check{\mathbf{x}}-{\mathbf{x}}^{*}|$ is not easy to estimate. The only certainty is that $|\check{\mathbf{x}}-{\mathbf{x}}^{*}|$ is related to the time step $\Delta t$. Similarly, the error of $|u\left({\mathbf{x}}^{*}\right)-u_h\left({\mathbf{x}}^{*}\right)|$ encompasses errors arising from the projection of mesh points onto the tangent plane ${\mathrm{T}}_c$ , in addition to interpolation errors. Consequently, the presence of $|\check{\mathbf{x}}-{\mathbf{x}}^{*}|$ and the projection error in $|u\left({\mathbf{x}}^{*}\right)-u_h\left({\mathbf{x}}^{*}\right)|$ results in a deterioration in the convergence order of the CFEM when compared to the MCFEM.

4. Numerical Examples

In this section, we evaluate the precision of the MCFEM under both diffusion-dominated and convection-dominated conditions. Then, the impact of variation of curvature and multi-connected surfaces on the MCFEM is verified. To simulate the phenomenon of pollutant injection on torus, a discontinuous source term problem is employed. Furthermore, the nonlinear convection velocity’s impact on the MCFEM is evaluated by applying the Burgers equation on a peanut-shaped surface. Finally, the impact of random initial conditions on our method is confirmed by employing the convection Allen-Cahn equation on other multi-connected surface. The $L^2$-errors (denoted by ${\mathrm{Err}}_{L^2}=\left|\right|u\left(T\right)-u_h^N{\left|\right|}_{L^2\left({\Gamma }_h\right)}$) and the $H^1$-errors (denoted by ${\mathrm{Err}}_{H^1}=\left|\right|u\left(T\right)-u_h^N{\left|\right|}_{H^1\left({\Gamma }_h\right)}$) of the numerical solutions are computed, respectively.

4.1. Accuracy Test on the Sphere

Initially, we will evaluate the spatial accuracy of the MCFEM and compare it with the CFEM, assuming a time step $\Delta t\approx h^2$. Consider the CRD equation (4) on a sphere

$$\psi (x,y,z)=x^2+y^2+z^2-0.25=0,$$

(41)

where the reaction coefficient $\mu$ is set to 1 and the convection velocity $\mathit{\beta }$ is set to $[0,0,0.5]$. The exact solution can be expressed as follows,

$$u(x,y,z,t)=t^2(1-tanh\left(\frac{z}{\sqrt{ϵ}}\right)).$$

(42)

When the diffusion parameter $ϵ\ll {\left|\right|\beta \left|\right|}_{L^2(\Gamma )}$, the exact solution u is discontinuous at the equator of sphere (42), resulting in a convection-dominated (singular perturbation) problem. To investigate this phenomenon, we set the diffusion parameter $ϵ$ to 1, 1E-1, 1E-2, 1E-3 and 1E-4 respectively. The corresponding results are presented in Table 1 to Table 5.

Table 1. The errors of different method with $ϵ=1$ and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$.

Table 1. The errors of different method with $ϵ=1$ and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$.

Table 2. The errors of different method with $ϵ=1$E-1 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-1.

Table 2. The errors of different method with $ϵ=1$E-1 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-1.

Table 3. The errors of different method with $ϵ=1$E-2 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-2.

Table 3. The errors of different method with $ϵ=1$E-2 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-2.

Table 4. The errors of different method with $ϵ=1$E-3 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-3.

Table 4. The errors of different method with $ϵ=1$E-3 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-3.

Table 5. The errors of different method with $ϵ=1$E-4 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-4.

Table 5. The errors of different method with $ϵ=1$E-4 and $\frac{2ϵ}{{\left|\right|\mathit{\beta }\left|\right|}_{L^2(\Gamma )}^2}=1.27$E-4.

The $L^2$-errors and $H^1$-errors of the MCFEM in Table 1 show the same trend as that of the CFEM. When $h>2.67$E-2, the mesh size h is not fine enough to ignore the existence of geometric errors. Consequently, the outcomes for $h\ge 2.67$E-2 deviate from anticipated results. As the mesh size h diminishes, the $L^2$-errors convergence order of the MCFEM attains 2, while the $H^1$-errors convergence order attains 1. This observation indicates that the MCFEM has second-order spatial accuracy when diffusion is dominant.

With the decrease of parameter $ϵ$, the diffusion of the equation begins to weaken and convection gradually dominates. The $L^2$-errors and $H^1$-errors of the MCFEM as shown in Table 2 and Table 3 are smaller than those of the CFEM. When $ϵ=1$E-3, the discontinuity of the exact solution u is obviously enhanced, resulting in a noticeable increase in the error of the MCFEM compared to $ϵ>1$E-3. Fortunately, the MCFEM still maintains second-order spatial accuracy in this case. Comparing with the MCFEM, the $L^2$-errors convergence order of the CFEM fails to reach 2 when $ϵ=1$E-1. This decrease occurs as a result of the discrepancy in u on the same normal vector when convection-dominated. The projection error in the CFEM’s reconstruction method cannot be disregarded. However, the MCFEM only involves Taylor expansion and surface finite element discretization without introducing other spatial errors. This is why our method ensures the $L^2$-errors convergence order is $O\left(h^2\right)$.

The continuity of the exact solution u is significantly compromised when the diffusion parameter $ϵ=1$E-4 as opposed to $ϵ=1$E-3. Consequently, the numerical solutions of the MCFEM and the CFEM under coarse mesh have obvious oscillation, as shown in Table 5. Analogous to the CFEM, the MCFEM demonstrates stability solely when $h=1.32$E-2. This observation underscores the influence of continuity on the mesh requirements for the MCFEM. Thus, it becomes imperative to employ a finer mesh size h to ensure the efficacy of the MCFEM when the diffusion parameter $ϵ$ is very small.

To visually reveal the distinctions between our proposed the MCFEM and the CFEM, we present a comparison of numerical solutions and error for various methods at $h=2.67$E-2, as illustrated in Figure 3 and Figure 4. As depicted in Figure 3, the $L^2$-errors of both the MCFEM and the CFEM exhibit an upward trend over time. Prior to $t=0.04$, the $L^2$-errors of the MCFEM is equivalent to that of the CFEM. With the increase of time, the $L^2$-errors growth rate of the CFEM is obviously faster than that of the MCFEM. This observation indicates that the accumulation of errors over time is significantly smoother for the MCFEM compared to the CFEM. The errors at the final moment of the MCFEM and the CFEM are shown in the subgraph $(c,e)$ in Figure 4. we can see that the error distribution of the MCFEM is sparsely concentrated near the equator and is an order of magnitude smaller than that of the CFEM. In contrast, the CFEM produces a narrow error band near the equator with slight oscillations above it. These findings suggest that the MCFEM produces less error than the CFEM once it reaches stability.

Furthermore, we also considered the curvature variation surfaces

$$\begin{array}{c}\hfill \mathrm{tooth}:\psi (x,y,z)=x^4+y^4+z^4-(x^2+y^2+z^2)=0,\end{array}$$

(43)

and multi-connected surface

$$\begin{array}{c}\hfill \mathrm{torus}:\psi (x,y,z)={(0.5-\sqrt{x^2+y^2})}^2+z^2-0.1^2=0.\end{array}$$

(44)

The exact solution will be modified by

$$u(x,y,z,t)=e^{-t}\frac{xy}{\pi }arctan\left(\frac{z}{\sqrt{ϵ}}\right),$$

(45)

while maintaining the convection velocity $\mathit{\beta }$ and reaction coefficient $\mu$ unchanged. The diffusion parameter $ϵ$ is set to 1E-3 and the time step $\Delta t$ is set to $2h^2$. Additionally, the exact solution of u at $T=0.5$ is simulated using the MCFEM and the CFEM on a tooth and a torus, respectively. The corresponding results are shown in Figure 5 and Figure 6.

As depicted in Figure 5 and Figure 6, the error is centered at $z=0$, aligning with anticipated results. Our proposed the MCFEM demonstrates its efficacy in handling surfaces with varying curvatures and multi-connected topologies. The error of both the MCFEM and the CFEM suggest that the stabilized the MCFEM outperforms the CFEM in terms of accuracy.

4.2. The Discontinuous Source Term Problem on Torus

Here, the MCFEM will be utilized to simulate the movement of pollutants on a torus, which is tantamount to resolving the convection-dominated problem that contains discontinuous source terms. The level set function expression of the torus is indicated in equation (44). In order to maintain stability within the MCFEM, it is imperative to restrict the mesh size to $h=1.25$E-2 and the time step to $\Delta t=1$E-3. Additionally, the reaction coefficient should be fixed at $\mu =1$, the parameter at $ϵ=1$E-3, and the convection-velocity at $\mathit{\beta }=[-y,x,0]$. Furthermore, the pollutant is yet to be introduced into the torus at time $t=0$, therefore an initial value of $u_0=0$ is selected. Four points ${\left\{{\mathbf{x}}_s\right\}}_{s=1}^4$ on the torus are arbitrarily selected, and the pollutants are continuously injected at these points ${\left\{{\mathbf{x}}_s\right\}}_{s=1}^4$ respectively to create a discontinuous source term

$$f=\left\{\begin{array}{cc}5,\hfill & |\mathbf{x}-{\mathbf{x}}_s|<0.03,s=1,2,3,4,\mathbf{x}\in \Gamma ,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(46)

According to the findings presented in Figure 7, pollutants form four stable regions on the torus under the influence of discontinuous source terms and convection velocity. This observation illustrates that the numerical solutions acquired through the MCFEM demonstrate comparable physical phenomena with the ones obtained from the CFEM in [32].

4.3. The Burgers Equation on Peanut-Shaped Surface

To investigate the impact of nonlinear problems on our methods, we selected the convection velocity $\mathit{\beta }=[u,0,0]$. Assigning the diffusion parameter $ϵ=1$E-3, reaction coefficient $\mu =0$, and setting the source term to $f=0$, the problem transforms into a typical Burgers problem on surfaces. Without loss of generality, a peanut-shaped surface

$$\psi (x,y,z)=({(2x-1)}^2+4y^2+4z^2)({(2x+1)}^2+4y^2+4z^2)-1.5=0,$$

(47)

is selected and the initial condition is

$$u_0=\left\{\begin{array}{cc}-sin\left(2\pi x\right),\hfill & \left|y\right|\le 0.5,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(48)

The corresponding mesh size $h=1$E-2 and time step $\Delta t=1$E-3. Since the convection velocity $\mathit{\beta }$ depends on time, the velocity $\mathit{\beta }$ at the current moment is approximated by the velocity $\mathit{\beta }$ at the previous moment.

The findings are presented in Figure 8. As time progresses, the numerical solution displays a marked modification at the centre of the peanut. To depict this trend visually, we projected the numerical solution $u_h^n$, where $\left|y\right|<4\mathrm{E}-4$, onto the X-axis, and the result is illustrated in Figure 9. The wave clearly propagates forward over time, and the gradient near $x=0$ exhibits progressive increments. The results of the calculations bear similarity to the one-dimensional case [35], implying the applicability of the MCFEM in addressing nonlinear convection-dominated problem on surfaces.

4.4. The Convection Allen-Cahn Equation on Multi-Connected Surface

In this example, we will analyze the impact of random initial conditions on the MCFEM by the convection Allen-Cahn equation

$${\partial }_{\tau }u-\epsilon {\Delta }_{\Gamma }u=f\left(u\right),$$

(49)

on a multi-connected surface

$$\psi (x,y,z)={(x^2+y^2-4)}^2+{(x^2+z^2-4)}^2+{(z^2+y^2-4)}^2+{(x^2-1)}^2+{(y^2-1)}^2+{(z^2-1)}^2-15=0.$$

(50)

The convection velocity $\mathit{\beta }=[0,0,-2]$, the diffusion parameter $\epsilon =1$E-3 and nonlinear function $f\left(u\right)=\frac{1}{\epsilon }(u^3-u)$ are selected. The initial condition is randomly chosen within the range of -0.1 to 0.1, as depicted in Figure 10. We define the discrete free energy

$$E_h^n={\int }_{{\Gamma }_h}\epsilon {\left|{\nabla }_{{\Gamma }_h}{u}_h^n\right|}^2+F\left(u_h^n\right)d\sigma ,$$

(51)

where $F\left(u_h^n\right)=\frac{1}{\epsilon }{({\left(u_h^n\right)}^2-1)}^2$. Additionally, the nonlinear function $f\left(u_h^n\right)$ can be approximated as $f\left(u_h^{n-1}\right)+\frac{2}{\epsilon }(u_h^n-u_h^{n-1})$.

The decrease of energy is a well-established properties of the convection Allen-Cahn equation. To observe the energy variation of the numerical solution, we control the mesh size $h=1.08$E-1 and the time step $\Delta t=1$E-2 to obtain Figure 11. The presented data in Figure 11 demonstrates a decrease in discrete energy over time, ultimately reaching a state of stability at $t=56$. This observation indicates that the random initial condition does not significantly impact the effectiveness of the MCFEM. To visually depict the progression from the initial condition to the steady state, numerical solutions at various time points within the range of $[0,100]$ were extracted and uniformly rotated. The outcomes are illustrated in Figure 12, revealing that the time trend of phase separation is consistent with the results observed in Figure 11.

5. Conclusions

This paper introduces a modified characteristic finite element method that exhibits second-order spatial accuracy. Our method employs Taylor expansion to reconstruct the solution beyond the surface in the characteristic direction. In contrast, the CFEM’s [32] reconstruction method introduces additional spatial errors, resulting in a lower spatial convergence order compared to our method. The reason for this phenomenon is that the definition domain of the solution has been extended from the surface to the neighborhood containing the surface by the characteristic finite element method. Consequently, the solutions along the same normal vector to be unequal by default, which is contrary to [26]. Despite the superior spatial accuracy of our proposed the MCFEM in comparison to the CFEM, this advantage comes at the expense of stability. The reason for this sacrifice in stability is that our reconstruction method transforming the characteristic finite element method to an explicit-implicit method.