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

Abstract

We present a modified characteristic finite element method that exhibits second-order spatial accuracy for solving convection–reaction–diffusion equations on surfaces. The temporal direction adopted the backward-Euler method, while the spatial direction employed the surface finite element method. In contrast to regular domains, it is observed that the point in the characteristic direction traverses the surface only once within a brief time. Thus, good approximation of the solution in the characteristic direction holds significant importance for the numerical scheme. In this regard, Taylor expansion is employed to reconstruct the solution beyond the surface in the characteristic direction. The stability of our scheme is then proved. A comparison is carried out with an existing characteristic finite element method based on face mesh. Numerical examples are provided to validate the effectiveness of our proposed method.

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, pattern formation [2,3], surfactant transportation [4], evolution of colonies on irregular surfaces [5], and cell migration [6,7]. 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 [8,9,10,11,12] and mesh-based methods [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]. 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 the mesh-based methods for solving PDEs on surfaces. The mesh generation include two prevalent strategies: embedding the surface in the narrow-band domain [16,17,18,19,20,21,22,23] and directly discretizing the surface [25,26,27,28,29,30].

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. [22], who have presented error estimates for a trace finite element method with streamline-upwind/Petrov–Galerkin stabilization. This represents the first residual-based stabilization method for convection–diffusion equations on surfaces. However, the accuracy and stability of this method depend on the careful selection of the stabilization 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 in their study, which is symmetric stability terms [26]. Recently, Xiao et al. developed gradient recovery-based adaptive techniques in [27]. Jin et al. [28] subsequently utilized their work [27] to address the convection-dominated problem on surfaces using mixed finite element methods and provided theoretical analysis.

The characteristic finite element method [23,29,31] is easier to operate than the above methods [22,26,27,28]. 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 [31]. Unlike the regular domains [31], 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 [23]. However, it presents challenges for the face mesh-based method in Section 4 of [25]. In response, Xiao et al. [29] introduced a characteristic finite element method (CFEM) that relies on the face mesh. By incorporating mass lumping, 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 method. 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 deterioratation 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 our findings 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 obtained with

$${\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

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

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

where $\mathbf{\nu }$ is the co-normal of Γ, and $dA$ and $d\sigma$ are the measures 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 with

$${\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 }$ are 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

$$\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,$$

(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,\dots ,N$, and 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 that there is a characteristic direction $\tau$ in $U(\Gamma ,\delta )\times [t^{n-1},t^n]$ and the point $\left(\mathit{\chi }\right(t),t)$ on direction $\tau$ satisfy

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

(7)

Integrating the characteristic Equation (7) over $[t^{n-1},t^n]$, the backtracking characteristic point obtained is as follows:

$$\mathit{\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(\mathit{\chi }\right(t),t)$ is

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

(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 analyze 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 [16,17,18,19,20,21,22,23] is used, the solution beyond the surface can be reconstructed using interpolation. However, the face mesh in Section 4 of [25] 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 used by Xiao et al. [29], 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 [29] 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 reconstructing 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 obtained with $\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 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:

$$\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.$$

(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 )}\hfill \\ \hfill & =\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 \\ \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 \\ \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 with

$$\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 )}\hfill \\ \hfill & =\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 completes the proof. □

Let $u^n$ be 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,\dots ,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 vs.\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,\dots ,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}{c}\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\}\\ \hfill \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}dt\right),\end{array}$$

(19)

for all n.

Proof. 

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

$${(\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 }.$$

(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\hfill \\ \hfill & =2{(u^{n-1}-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^{n-1}+\Delta t{f}^n,u^n)}_{\Gamma }\hfill \\ \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 \\ \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 }+\Delta t^2\left|\right|f^n{\left|\right|}_{L^2(\Gamma )}^2+\left|\right|u^n{\left|\right|}_{L^2(\Gamma )}^2\hfill \\ \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\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}$$

(23)

It follows (5) that the bound

$$-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 },$$

(24)

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

$$\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.$$

(25)

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

$$\begin{array}{cc}\hfill & \left|\right|u^n-\Delta t{\mathit{\beta }}_{\Gamma }·{\nabla }_{\Gamma }{u}^n{\left|\right|}_{L^2(\Gamma )}^2\hfill \\ \hfill & ={(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 of 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\hfill \\ \hfill & \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 \\ \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 holds:

$$\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),$$

(28)

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 is composed of plane open triangles $K_j$ with edge $\partial K_j$ and vertexes ${\mathbf{x}}_l,l=1,\dots ,N_h^v$. The point ${\mathbf{x}}_l$ is also the vertex points of the 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,\dots ,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 [25], 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 obtained with $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 obtain the MCFEM using Equation (10): For $n=1,2,\dots ,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 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$ and $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|[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)\right|\le \kappa \forall h<h_{\kappa },$$

where $[[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h},{\Gamma }_h}]]\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 & [[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)\hfill \\ \hfill & =\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 & [[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)\hfill \\ \hfill & =\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 & [[\mathbf{\nu }·{\mathit{\beta }}_{{\mathcal{P}}_{{\Gamma }_h,{\Gamma }_h}}]]\left({\mathcal{P}}_{{\Gamma }_h}\left({\mathbf{x}}_0\right)\right)\hfill \\ \hfill & =\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}{c}\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\}\\ \hfill \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+{\int }_0^T\left|\right|f_{{\mathcal{P}}_{{\Gamma }_h}}\left(t\right){\left|\right|}_{L^2\left({\Gamma }_h\right)}^{2}dt\right),\end{array}$$

(33)

for $n=1,2,\dots ,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\hfill \\ \hfill & =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 \\ \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}\hfill \\ & \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+{\int }_{t^{n-1}}^{t^n}\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}\hfill \\ \hfill & \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}\hfill \\ \hfill & ={\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 \\ \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 \\ \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$, and ${\mathbf{\nu }}_{\partial K_j}$ is the 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 measures 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 that $K_j$ is a triangular plane and 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 with

$$\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}\hfill \\ \hfill & =\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 \\ \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 \\ \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 obtain

$$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.$$

(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 Section 3.1 is different from that employed in the CFEM [29], 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 Section 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 neighborhood $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 surface finite element method described in [25].

Employing the notations given in Section 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 in 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 of $\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 (\frac{z}{\sqrt{ϵ}})).$$

(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, 1 $\times {10}^{-1}$, 1 $\times {10}^{-2}$, 1 $\times {10}^{-3}$ and 1 $\times {10}^{-4}$, respectively. The corresponding results are presented in

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

h	MCFEM				CFEM
	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate
2.04 $\times {10}^{-1}$	2.42 $\times {10}^{-2}$	–	4.33 $\times {10}^{-2}$	–	2.30 $\times {10}^{-2}$	–	4.23 $\times {10}^{-2}$	–
1.08 $\times {10}^{-1}$	1.44 $\times {10}^{-2}$	0.82	2.70 $\times {10}^{-2}$	0.74	1.49 $\times {10}^{-2}$	0.69	2.80 $\times {10}^{-2}$	0.65
5.38 $\times {10}^{-2}$	3.63 $\times {10}^{-3}$	1.99	1.05 $\times {10}^{-2}$	1.37	3.85 $\times {10}^{-3}$	1.96	1.08 $\times {10}^{-2}$	1.37
2.67 $\times {10}^{-2}$	4.27 $\times {10}^{-4}$	3.05	4.57 $\times {10}^{-3}$	1.18	3.98 $\times {10}^{-4}$	3.23	4.57 $\times {10}^{-3}$	1.23
1.32 $\times {10}^{-2}$	1.05 $\times {10}^{-4}$	2.00	2.25 $\times {10}^{-3}$	1.01	1.09 $\times {10}^{-4}$	1.85	2.26 $\times {10}^{-3}$	1.00

,

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

h	MCFEM				CFEM
	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate
2.04 $\times {10}^{-1}$	2.65 $\times {10}^{-2}$	–	1.34 $\times {10}^{-1}$	–	2.19 $\times {10}^{-2}$	–	1.31 $\times {10}^{-1}$	–
1.08 $\times {10}^{-1}$	1.68 $\times {10}^{-2}$	0.71	6.86 $\times {10}^{-2}$	1.05	2.01 $\times {10}^{-2}$	0.14	7.76 $\times {10}^{-2}$	0.82
5.38 $\times {10}^{-2}$	4.24 $\times {10}^{-3}$	1.99	3.19 $\times {10}^{-2}$	1.11	6.01 $\times {10}^{-3}$	1.74	3.58 $\times {10}^{-2}$	1.12
2.67 $\times {10}^{-2}$	4.75 $\times {10}^{-4}$	3.12	1.54 $\times {10}^{-2}$	1.04	1.04 $\times {10}^{-3}$	2.51	1.63 $\times {10}^{-2}$	1.12
1.32 $\times {10}^{-2}$	1.17 $\times {10}^{-4}$	2.00	7.59 $\times {10}^{-3}$	1.01	5.55 $\times {10}^{-4}$	0.89	8.15 $\times {10}^{-3}$	0.99

,

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

h	MCFEM				CFEM
	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate
2.04 $\times {10}^{-1}$	4.15 $\times {10}^{-2}$	–	7.20 $\times {10}^{-1}$	–	4.77 $\times {10}^{-2}$	–	7.31 $\times {10}^{-1}$	–
1.08 $\times {10}^{-1}$	1.90 $\times {10}^{-2}$	1.22	3.39 $\times {10}^{-1}$	1.18	3.48 $\times {10}^{-2}$	0.49	4.25 $\times {10}^{-1}$	0.85
5.38 $\times {10}^{-2}$	4.73 $\times {10}^{-3}$	2.00	1.64 $\times {10}^{-1}$	1.05	1.54 $\times {10}^{-2}$	1.18	2.29 $\times {10}^{-1}$	0.89
2.67 $\times {10}^{-2}$	6.64 $\times {10}^{-4}$	2.80	8.04 $\times {10}^{-2}$	1.01	6.59 $\times {10}^{-3}$	1.21	1.17 $\times {10}^{-1}$	0.96
1.32 $\times {10}^{-2}$	1.63 $\times {10}^{-4}$	2.00	3.94 $\times {10}^{-2}$	1.02	3.43 $\times {10}^{-3}$	0.93	5.99 $\times {10}^{-2}$	0.95

,

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

h	MCFEM				CFEM
	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate
2.04 $\times {10}^{-1}$	1.03 $\times {10}^{-1}$	–	2.86 $\times {10}^0$	–	1.09 $\times {10}^{-1}$	–	2.65 $\times {10}^0$	–
1.08 $\times {10}^{-1}$	4.95 $\times {10}^{-2}$	1.15	2.40 $\times {10}^0$	0.27	7.23 $\times {10}^{-2}$	0.65	2.15 $\times {10}^0$	0.33
5.38 $\times {10}^{-2}$	1.18 $\times {10}^{-2}$	2.07	1.32 $\times {10}^0$	0.92	4.10 $\times {10}^{-2}$	0.82	1.55 $\times {10}^0$	0.47
2.67 $\times {10}^{-2}$	1.47 $\times {10}^{-3}$	2.97	4.54 $\times {10}^{-1}$	1.50	2.14 $\times {10}^{-2}$	0.92	8.29 $\times {10}^{-1}$	0.89
1.32 $\times {10}^{-2}$	3.37 $\times {10}^{-4}$	2.10	2.21 $\times {10}^{-1}$	1.03	1.12 $\times {10}^{-2}$	0.93	4.31 $\times {10}^{-1}$	0.93

and

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

h	MCFEM				CFEM
	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate	${\mathbf{Err}}_{{\mathit{L}}^{\mathbf{2}}}$	Rate	${\mathbf{Err}}_{{\mathit{H}}^{\mathbf{1}}}$	Rate
2.04 $\times {10}^{-1}$	1.04 $\times {10}^{-1}$	–	5.90	–	1.54 $\times {10}^{-1}$	–	5.33	–
1.08 $\times {10}^{-1}$	1.53 $\times {10}^{-1}$	−0.61	5.49	0.11	1.07 $\times {10}^{-1}$	0.56	5.06	0.80
5.38 $\times {10}^{-2}$	8.20 $\times {10}^{-2}$	0.90	8.32	−0.60	7.84 $\times {10}^{-2}$	0.45	5.41	−0.97
2.67 $\times {10}^{-2}$	2.75 $\times {10}^{-2}$	1.56	5.94	0.48	4.64 $\times {10}^{-2}$	0.75	4.66	0.22
1.32 $\times {10}^{-2}$	2.78 $\times {10}^{-3}$	3.27	1.57	1.90	2.41 $\times {10}^{-2}$	0.94	2.65	0.81

.

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$ $\times {10}^{-2}$, the mesh size h is not fine enough to ignore the existence of geometric errors. Consequently, the outcomes for $h\ge 2.67$ $\times {10}^{-2}$ deviate from the 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 in 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$ $\times {10}^{-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$ $\times {10}^{-3}$. Fortunately, the MCFEM still maintains second-order spatial accuracy in this case. Compared with the MCFEM, the $L^2$ errors convergence order of the CFEM fails to reach 2 when $ϵ=1$ $\times {10}^{-1}$. This decrease occurs as a result of the discrepancy in u on the same normal vector when it is 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 that 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$ $\times {10}^{-4}$ as opposed to $ϵ=1$ $\times {10}^{-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$ $\times {10}^{-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 MCFEM and CFEM, we present a comparison of numerical solutions and error for various methods at $h=2.67$ $\times {10}^{-2}$, as illustrated in Figure 3 and Figure 4. As depicted in Figure 3, the $L^2$ errors of both the MCFEM and 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 in 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 with

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

(45)

while maintaining the convection velocity $\mathit{\beta }$ and keeping the reaction coefficient $\mu$ unchanged. The diffusion parameter $ϵ$ is set to 1 $\times {10}^{-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 the anticipated results. Our proposed MCFEM demonstrates its efficacy in handling surfaces with varying curvatures and multi-connected topologies. The errors of both the MCFEM and the CFEM suggest that the stabilized 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$ $\times {10}^{-2}$ and the time step to $\Delta t=1$ $\times {10}^{-3}$. Additionally, the reaction coefficient should be fixed at $\mu =1$, the parameter at $ϵ=1$ $\times {10}^{-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 are 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 [29].

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$ $\times {10}^{-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 as follows:

$$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$ $\times {10}^{-2}$ and time step $\Delta t=1$ $\times {10}^{-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\times {10}^{-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 [32], 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 using 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$ $\times {10}^{-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 as

$$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 in energy is a well-established property of the convection Allen–Cahn equation. To observe the energy variation in the numerical solution, we control the mesh size $h=1.08$ $\times {10}^{-1}$ and the time step $\Delta t=1$ $\times {10}^{-2}$ to obtain Figure 11. The presented data in Figure 11 demonstrate 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 [29] 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 with the characteristic finite element method. Consequently, the solutions along the same normal vector remain unequal by default, which is contrary to [25]. Despite the superior spatial accuracy of our proposed 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 transforms the characteristic finite element method to an explicit–implicit method.