A fast method for solving time-dependent nonlinear convection diffusion problems

1.   Introduction

In this paper, we introduce a fast IMEX (implicit-explicit) scheme for the 2D unsteady nonlinear convection-diffusion problems. Our model problem is as follows

with the Dirichlet boundary condition

and the initial value condition

where $\Omega\subset R^2$ is a bounded polygonal domain with Lipschitz continuity, $P(u)$ and $Q(u)$ are the convective field functions, $\varepsilon > 0$ is a given diffusion coefficient, and $T$ is the final time.

Note that we can rewrite the system of Eqs (1.1)–(1.3) into the vector form as

with the Dirichlet boundary and initial value conditions

where $F(u) = (P(u), Q(u))^T$, ${\bf{x}} = (x, y)^T$, ${u} = u({{\bf x}}, t)$.

Generally, some physical phenomenon, where the energy or other physical substances are transferred in physical systems as a result of diffusion and convection processes, can be described by nonlinear convection-diffusion Eqs (1.1)–(1.3). These systems are of reasonably practical significance since they can describe varieties of physical phenomenon, which appear in a plenty of applications and symbolize mathematical models for a great deal of physical processes in mechanics of fluids, astrophysics, groundwater flow, meteorology, semiconductors and reactive flows. In particular, for the certain case when $f = 0$, (1.1)–(1.3) are scalar convection-diffusion equations, which usually arise in various applications, with the range including from turbulence model ^[1,2], through traffic flow ^[3], to two-phase flow in a porous medium ^[4].

For the general nonlinear convection-diffusion Eqs (1.1)–(1.3), a group of high order finite difference numerical schemes were raised based on Euler implicit scheme by Cecchi and Pirozzi ^[5], a discontinuous Galerkin finite element discretization was introduced to study the effect of numerical integration by Sobotikova et al. ^[6], and a Lattice Boltzmann model was shown by Shi and Guo ^[7]. For the certain case $F = 0$, namely the nonlinear conservative laws, Karlsen et al. ^[8] used the idea of operator splitting to decompose the convection part and the diffusion part of the problem, and provided a dimensional splitting method to simplify the m-dimensional convection equation to a series of one-dimensional equations. After that, Karlsen et al. ^[9] presented a corrected operator splitting method to deal with much more complicated problems with general flux functions and initial value. Many central schemes were also constructed to achieve high resolution of the numerical solutions; See for instance Nessyahu and Tadmor ^[10], Jin and Xin ^[11], Kurganov and Tadmor ^[12]. Recently, some high order Maximal principle-Preserving finite difference Rungeutta WENO (weighted essentially nonoscillatory) or DG (discontinuous Galerkin) schemes were constructed; See Jiang and Xu ^[13], and Xiong et al. ^[14] and the references therein.

As explained by Chertock and Kurganov ^[15], such problems are parabolic and that's why even for discontinuous initial value, they can be admitted global smooth solutions. In this paper, we combine the continuous finite element scheme and operator splitting approach to propose a fast method for solving these problems. It is obvious that the treatment of nonlinear convection term in partial differential equations like (1.1)–(1.3) is one of the key procedures, which affects the efficiency of the numerical methods. An important class of numerical methods about this issue is referred as the IMEX method; See for instance Akrivis ^[16] for general nonlinear parabolic problems, Long and Chen ^[17] for nonlinear convection diffusion equations, as well as Zhang et al. ^[18] for the slight simple non-stationary thermal convection problems.

In this paper, first of all, at each time step, the initial coupled problem is decomposed into convection and diffusion subproblems by employing the same idea of operator splitting with ^[8] (see also ^[19] for much general parabolic problems), and much importantly, make use of a specific Taylor expansion to deal with the nonlinear term in convection subproblems, and eventually provide an IMEX scheme for solving the nonlinear convection diffusion problems. Besides, for the sake of relaxing the possible poor stability condition that results from the explicit operation of the convection problem, a multistep technique is proposed. There are some major superiorities of the present method. Firstly, at each time step, since the self-adjointness and coerciveness of diffusion subproblems, there are a lot of forceful preconditioned iterative solvers, like preconditioned conjugate gradient solver, which can be used to solve them efficiently. Secondly, it is easy for us to utilize the mass lumping method to explicitly solve the mass matrices for spatial linear elements because of the linearity of convection subsystems. When we use the finite element method to discretize the equations, all the major stiffness matrices keep invariable in the time stepping process, that is the reason why the present method is notably fast for solving the unsteady nonlinear systems.

The paper is formed as follows. We propose the idea of construction of the present semi-discrete IMEX method via operator splitting, which is shown in Section 2. After that, by means of the spatial finite element discretization, a further representation of the algorithm of full discretization is illustrated. Meanwhile, the superiorities of the proposed method are illuminated by some remarks as well as the multistep scheme. As the numerical study, several numerical experiments on Burgers type and Buckley-Leverett type equations and discussions are presented in Section 3.

2.   The IMEX scheme

In this section, we propose the semi-discrete and its corresponding fully discrete IMEX algorithms. Without loss of generality, we uniformly split the time interval to make temporal discretization. Fix positive integer $N\ge1$, set $0 = {t^0} < {t^1} \cdots < {t^{N-1}} < {t^N} = T$, with ${t^n} = n \Delta t$, ${t^{n+\frac12}} = {t^n}+ \Delta t/2$, and the step size $\Delta t = T/N$.

Firstly, the nonlinear convection term $\nabla\cdot(F({u}))$ and the linear diffusion term $\varepsilon\Delta u$ are separated. We can write it in the form of operators described below

Where, $A_1u$ is the nonlinear part of the convection diffusion equations. In order to compute the Eqs (1.4) and (1.5), we can directly utilize numerous operator splitting methods, see, e.g., ^[8] by Karlsen et al. and [19, Chapter Ⅱ] by R. Glowinski and among the references. In this paper, we employ the following time splitting technique

on the time interval $[t^n, t^{n+1}]$ to solve the original nonlinear system (1.4)–(1.5). When it comes to the time splitting, the certain Euler central difference and backward schemes are utilized to discretize (2.1) and (2.2), seperately,

where $F(v({{\bf x}}, t^{{n+\frac12}}))$ is approximated by $F^{n+\frac12}$, and the solution $u({{\bf x}}, t)$ to the system (1.4)–(1.5) at time instances $t = t^n$ and $t^{n+1}$ are approximated by $u^n$ and $u^{n+1}$, correspondingly.

It worth noting that the performance of the resulting algorithm, such as the stability, convergence and efficiency, would be significantly affected by the numerical treatment of $F^{n+\frac12}$ in the iterative scheme (2.3)–(2.4). For instance, seting $F^{n+\frac12} = (F(v^{n+1})+F({v}^{n}))/2$ in (2.3) results in a Crank-Nicolson scheme, and generally speaking, we still get a nonlinear resulting system. If we take $F^{n+\frac12}$ as an explicit linearization, i.e., $F^{n+\frac12} = F(v^n) = F(u^n)$, Eq (2.3) will generate a Forward Euler scheme with strong instability generally. Without choosing $F^{n+\frac12}$ as these two ways, we linearize $F^{n+\frac12}$ at the time instance $t = t^n$ by proposing a new unique Taylor expansion method. More accurately, as a result of the fact that Eq (2.1) is discretized into (2.3), then we can access through Taylor expansion:

here we use the definition $F_u = \frac{d F(u)}{d u}$. It is noted that throughout this paper, we assume that $F(u)\in C^1(R)$ in order to construct our present numerical algorithm; See Remark 1 in Reference ^[6] for more details.

In term of (2.1) and (2.3) we obtain that

and

then a linearization of $F^{n+\frac12}$ is derived as

Denoting $\eta^n = F({u}^n)-\frac12{\Delta t} F_u(u^n)\nabla\cdot (F(u^n))$, the following time splitting scheme is the newly proposed one that we now obtain:

Meanwhile, the boundary condition for these Eqs (2.5) and (2.6) is need to be fulfilled. When solving $u^{n+1}$ in the subsystem (2.6), further a temporal approximation of diffusion problem (2.2), we can naturally inherit a Dirichlet boundary condition on $\Gamma = \partial\Omega$ for ${u^{n+1}}$ from (1.5), namely

Since $u_*^{n+1}$ in the subsystem (2.5) is a nontrivial temporal approximation of convective Eq (2.1), we should only impose partial boundary condition. Here, a similar artificial inflow boundary condition is introduced

where $\Gamma_t^-$ is the inflow boundary named by

and the outward normal vector to $\Gamma$ is expressed by ${{\bf n}}$.

For simplicity in the rest of the paper, let us denote by $u_b^{n+1} = u_b(\cdot, t^{n+1})$, and functional space is defined as $X = { H}^1({\Omega})$. Besides, terms (2.5) and (2.6) are respectively multiplied by a testing functions $v\in X$ with $v|_{\Gamma_{n+1}^-} = 0$ and ${v}|_{\Gamma} = 0$, (as the options of boundary conditions (2.7) and (2.8)) correspondingly, and integrated by parts to obtain the variational formulations. Finally, we can describe the time splitting method as the following:

Temporal Semi-approximation Scheme:

$\bf{1.}$ Initialize value: Let ${u}^0 = {u}_0$. $\bf{2.}$ For any $n\ge 0$, supposing ${u}^n$ is already known, find ${u}_{*}^{n+1}$ and ${u}^{n+1}$ to solve

Convection subsystem: Compute ${u}_{*}^{n+1}\in X$, ${u}_{*}^{n+1}|_{\Gamma_{n+1}^-} = u_b^{n+1}$, solve

for any ${v} \in X$, with ${v} = 0$ on ${\Gamma_{n+1}^-}$;

and

Diffusion subsystem: Compute ${ u}^{n+1}\in X$, ${u}^{n+1}|_{\Gamma} = {u}_b^{n+1}$, such that

for any ${v} \in X$, with ${v} = 0$ on ${\Gamma}$.

$\bf{3.}$ Let $n+ = 1$, until $n > N$, else retrun to Step 2.

Next, the finite element discretization of the above time splitting method will be discussed with a shape-regular triangulation of ${\Omega}$, which is denoted by $\mathcal{T}_{h}$ and has mesh scale $h = \max_{K\in\mathcal{T}_{h}} {\rm{diam}}(K).$ Denote the finite element space with order $k\ge1$ as

that is a subspace of $X$ used to approximate the velocity. In the next segment, we simply offer a detailed discussion under the certain circumstance $k = 1$, which is the continuous piece-wise linear functions over $\mathcal{T}_{h}$. Apart from that, the standard nodal interpolation operator denoted by ${I}_h:X\rightarrow X_h$ is needed.

Next, the following scheme is our full (space-time) approximation approach.

Algorithm A: Space-time full approximation Scheme

$\bf{1.}$ Initialize value: Let ${u}_h^0 = { I}_h{u}_0$.

$\bf{2.}$ For any $n\ge 0$, supposing ${u}_h^n$ is already known, find ${u}_{h, *}^{n+1}$ and ${u}_h^{n+1}$ by solving

Convection subsystem: Denote $\eta_h^n = F(u_h^n)-\frac12F_u(u_h^n)\nabla\cdot(F(u_h^n))$; Compute ${u}_{h, *}^{n+1}\in X_h$, with ${u}_{h, *}^{n+1}|_{\Gamma_{n+1}^-} = { I}_hu_b^{n+1}$, solve

for any ${v}_h \in X_h$, with ${v}_h = 0$ on ${\Gamma_{n+1}^-}$;

and

Diffusion subsystem: Find ${ u}_h^{n+1}\in X_h$ with ${u}_h^{n+1}|_{\Gamma} = {I}_h{u}_b^{n+1}$, such that

for any ${v}_h \in X_h$, with ${ v}_h = 0$ on ${\Gamma}$. $\bf{3.}$ Let $n+ = 1$, until $n > N$, else retrun to Step 2.

Remark 1. For any given time step $\Delta t$, and for all $n\le N$, we can get some positive definite and permanently invariant stiffness matrices for (2.12).

Remark 2. Since the stiffness matrices of (2.11) are mass matrices as well as maintain invariable for fixed $\Gamma_n^- = \Gamma^-$ ($n\le N$), the subsystem can be explicitly solved. And particularly, when using linear finite element space $(k = 1)$, the mass-lumping technique works efficiently (see details in reference ^[20]).

Remark 3. Compared with the majority of existing techniques for computing (1.4)–(1.5), the distinct operation of nonlinear convection term would produce some unique superiorities of the present approach, as depicted in Remark 2. Nevertheless, as smaller time step size is required to realize convergence in this explicit linearization, the practical applications of our present approach would slightly decline. As a matter of fact, it is a good way to apply the multistep pattern of the approach to remedy this deficiency. Practically, when we deal with each diffusion subsystem (2.12) using time step size $\Delta t$, we can utilize a smaller $\delta t = \Delta t/m$, to complete the convection subsystem (2.11) for $m$ times.

More accurately, the multistep approach is summarized as follows.

Algorithm B: Multistep Scheme

$\bf{1.}$ Initialize value: Let ${u}_h^0 = { I}_h{u}_0$.

$\bf{2.}$ For any $n\ge 0$, supposing ${u}_h^n$ is already known, find ${u}_{h, *}^{n+1}$ and ${u}_h^{n+1}$ by solving

Convection subsystem: Let $u_{h, *}^n = u_h^n$; For any $i = 1, 2, \cdots, m$, solve

and compute ${u}_{h, *}^{n+\frac{i}{m}}\in X_h$, with ${u}_{h, *}^{n+\frac{i}{m}}|_{\Gamma_{{n+{i}/{m}}}^-} = { I}_hu_b^{{n+\frac{i}{m}}}$, such that

for any ${v}_h \in X_h$, with ${v}_h = 0$ on ${\Gamma_{n+i/m}^-}$;

and

Diffusion subsystem: Find ${ u}_h^{n+1}\in X_h$ with ${u}_h^{n+1}|_{\Gamma} = {I}_h{u}_b^{n+1}$, and compute

for any ${v}_h \in X_h$, with ${ v}_h = 0$ on ${\Gamma}$.

$\bf{3.}$ Let $n = n+1$, until $n > N$, else return to Step 2.

Remark 4. In this multistep algorithm, the inflow boundary $\Gamma_{n+\frac{i}{m}}^-$ for subsystem (2.13) is modified as $\Gamma_{n+\frac{i}{m}}^- = \left\{{{\bf x}}\in\Gamma:\, F_u(u_b({{\bf x}}, t^{{n+\frac{i}{m}}}))\cdot{{\bf n}}({{\bf x}}) < 0\right\},$ which is consistent with the previous definition $\Gamma_{n}^-$.

Remark 5. Actually, our multistep method is similar with the one used in wave-like equation strategy to compute convection advection equations (see e.g., Section 31.5 in ^[19] and the references therein). Likely, considering solving convection-dominated diffusion equations, the same multistep method was adopted and generally carried out, see for instance, References ^[21,22]. Moreover, we note that a multirate iterative scheme was proposed for solving poroelasticity model in a similar way of dealing with the Stokes and diffusion systems ^[23], more recently the multirate finite steps methods were provided to deal with the phenomena of different scales ^[24].

3.   Numerical experiments

This section aims to test the performance of the present scheme, including validness, stability and convergence. We use finite element package Freefem++ ^[25] to implement our scheme and carry out numerical experiments in the following.

3.1.   Burgers-type equations

This subsection focuses on the particular case $F(u) = (P(u), Q(u))^T = (\frac12 u^2, \frac12 u^2)^T$, and $f = 0$. Actually, these choices give rise to the equations of Burgers type.

Problem 1. The first example is a smooth problem on the unit square domain $\Omega = [0, 1]^2$, and the exact solution is given by $u(x, y, t) = 16x(1-x)y(1-y)(e^{10t}-1)/(e^{10}-1).$ The original and boundary conditions depend on (1.1)–(1.3) and the analytical solution. We choose $T = 1$, $\varepsilon = 10^{-8}$.

For this example, we can prove easily that for any $(x, y)\in \partial \Omega$ and $t\ge 0$, $F_u(u(x, y, t)) = (u, u)^T = (0, 0)^T$, as well as $F_u(u(x, y, t))\cdot{{\bf n}} = 0$, so the inflow boundary $\Gamma_n^{-}$ is always empty for the convective subsystem (2.11) at each time step $t = t^n$ for any positive integer $n$ satisfied the condition $t^n\le T$.

Firstly, we test the order of convergence with fixing time step $\Delta t = 1/2^{16}$ and varying the mesh scale from $1/4$ to $1/64$ step by step. In Table 1, the computational results are presented and we can observe an optimal convergence of rate 2 for the numerical solutions under $L^2$-norm error.

Next, we settle the spatial mesh with $h = 1/128$ and along with the time step size varying from $0.1$ to $0.1/2^6$. As it can be seen from Table 2, the Algorithm A (Single step scheme) can always derive the convergent results up to terminal time $T$. According to this table, we can clearly obtain an optimal convergence of rate 1 under the $L^2$-norm of the numerical solutions with regard to the time step size.

Problem 2. The second example is defined on the square $\Omega = [-1, 1]^2$, with the exact solution $u(x, y, t) = -(1-x^2)^2(1-y^2)^2\arctan\left(\frac{x+y+2-2t}{\nu}\right).$ The initial and boundary conditions depend on (1.1)–(1.3) and the analytic solution. Besides, the parameters are selected as $\varepsilon = 0.002$, $T = 1$ and $\nu = 0.5$.

Similarly with previous problem, for the Problem 2 we can check that, the inflow boundary $\Gamma_n^-$ for the convective subsystem (2.11) is awayls null.

In order to check the convergence order in spatial scale, a smaller fixed time step $\Delta t = 1/2^{16}$ is utilized for the present algorithm. As to this time step, Algorithm A can constantly get convergent solutions up to $T = 1.0$, and obtain an optimal convergence of rate 2 with respect to mesh scales by observing the numerical simulations under the $L^2$-norm error, which are shown in Table 3.

Now, we show numerical results in Table 4 with fixed spatial mesh $(h = 1/64)$, time step size from $0.1$ to $0.1/2^5$, and the same computational parameters. It displays the performance of the Algorithm B. Our observations for this testing case focus on the following two-folds: 1) An optimal convergence of order 1 under the convergence region with respect to the time step size is examined for the $L^2$-norm of the numerical simulations in time; 2) Algorithm A can get convergent solution for relatively large time step, for instance $\Delta t = 0.05$. Noting that for time step $\Delta t = 0.025$, the computational solution by Algorithm A shows apparently oscillation near the top region, which can be examined in the left contour plot of Figure 1. Compared with that, when using Algorithm B with multistep steps of $m = 4$ for convection subproblem, the smooth solution can be derived. This phenomena is observed in the right contour plot of Figure 1.

Next, we apply the Algorithm B for different time step sizes to compute the numerical solutions with fixing the multistep size $m = 64$, which are displayed in Table 5. We can notice that Algorithm B always get convergent and smooth solutions for reasonable large multisteps, namely $m = 64$ in this testing case. In sum, this multistep scheme behaves like an unconditional stable scheme in this aspect. Apart from that, we also observe the optimal convergence order with regard to time step size.

Moreover, at $\Delta t = 0.001$, the profiles of the solutions are shown in Figure 2, from which the smooth evolution of the numerical solutions at different time can be observed. Here, only the 3D plot of the numerical solutions at initial time and $t = 0.6, 0.7, 0.8, 0.9$ and terminal time $t = 1.0$ are shown.

3.2.   Buckley-leverett type equation

In this subsection, we will engage in the numerical simulations for two-dimensional Buckley-Leverett problem, which is a conservation system used to simulate two-phase flow in porous medium. The Buckley-Leverett problem or the Buckley-Leverett displacement describes an immiscible displacement process, such as the displacement of oil by water.

Problem 3. Consider the flux function in the form of $P(u) = u+(u-0.25)^3, Q(u) = -(u+u^2).$ The computational domain is $\Omega = [-0.5, 0.5]^2$, and the exact solution is $u(x, y, t) = 1/4+xy.$ The source is $f(x, y) = y+3x^2y^3-x(1.5+2xy)$, and the parameters are selected as $\varepsilon = 0.1, 0.01$, $T = 1$.

In this problem, we choose the same flux $F(u) = (P(u), Q(u))$ as the one from Example 2 of Reference ^[8]. For this flux, an exact solution with non-homogeneous boundary condition is constructed. Then, from these settings, since the exact solution is independent of time $t$, the present method can be interpreted as a pseudo-time method for solving nonlinear elliptic problems. In the present method (Algorithm A or B), we can easily check that $F_u\cdot{{\bf n}} < 0$ for the right and top boundaries, namely, $\{(x, y): x = 0.5, -0.5\le y\le 0.5; or -0.5\le x\le 0.5, y = 0.5\}$.

We apply Algorithm A to capture the numerical solutions up to $T = 1.0$ for two different viscosity parameters $\varepsilon = 0.1$ and $0.01$ with fixing a much smaller time step size $\Delta t = 1/2^{16}$. As it is shown in Table 6, the convergence results of the second and fourth columns correspond to the viscosity $\varepsilon = 0.1$ and $\varepsilon = 0.01$, respectively. From this table, an optimal convergence of rate 2 with respect to mesh scales is derived for the numerical simulations under $L^2$-norm for both viscosity parameters.

Next, the spatial mesh with $h = 1/128$ is fixed, and the time step size $0.1/2^k$ ($k = 1, 2, \cdots, 8$) are selected to check the performance of the Algorithm A or B. Firstly, for the case $\varepsilon = 0.1$, Table 7 shows that the single step scheme (Algorithm A) can always get the convergence results. From this table, the optimal first order convergence in time is verified under $L^2$-norm in time. Secondly, for the case $\varepsilon = 0.01$, as shown in the second column of Table 8, Algorithm A is convergent for $k\le 4$, and Algorithm B is then used by doubly increasing the multistep index $m$ each time until deriving the convergence numerical solution at the terminal time for each $k > 4$. Then for these selected time step sizes, we apply Algorithm B with fixed multistep index $m = 32$ to capture the numerical solutions, of which the approximate errors under $L^2$-norm at terminal time are listed in the three column of Table 8. We can give some observation based on this testing case: 1) An optimal convergence of rate 1 with regard to the time step size is derived for the numerical simulations under $L^2$-norm in time for Algorithm B with fixed index $m$; 2) For fixed time step size, Algorithm B with different multistep indices (noting that Algorithm B with $m = 1$ reduces to Algorithm A) can get very close approximate accuracy, which shows that the multistep scheme for convection subproblem in Algorithm B mainly plays a significant role in improving the stability of the algorithm, not the computational accuracy.

Problem 4. Consider the same flux function as in Problem 3 with the computed domain selected as $\Omega = [-2, 5]^2$. The initial value is chosen as $u(x, y, 0) = 1$ in $(x-0.25)^2+(y-2.25)^2 < 0.5$ and otherwise $u(x, y, 0) = 0$. The boundary condition is $u(x, y, t) = 0$, for any $(x, y)\in \partial \Omega$ and the right source is $f = 0$. The parameters are selected as $\varepsilon = 0.01$, $T = 1.0$.

This problem with discontinuous initial data was first given in ^[8] by Karlsen and Risebro, where only the convergence rates were plotted for different CFL numbers. As to this problem, similarly with the previous problem, we can deduce that the inflow boundaries for Eq (2.5) with condition $F_u\cdot{{\bf n}} < 0$ are the right and top boundaries, namely, $\{(x, y): x = -2, -2\le y\le 5; or -2\le x\le 5, y = 5\}$.

We execute Algorithm A under the discretization parameters $\Delta t = 0.001$ and $h = 1/64$ up to the terminal time $T = 1.0$. The contour plot of the numerical solution at $T = 1.0$ is displayed in Figure 3, where the contour values with Labels 1–18 correspond to the lines located from outside to inside in order. Besides, we displayed the three dimensional plot of the computational numerical approximations at different times in Figure 4. From this figure, we can observe the smooth evolution of the numerical solutions by Algorithm A, which verifies the stability of the proposed algorithm.

Problem 5. Consider the flux function of the form

The computed domain is $\Omega = [-1.5, 1.5]^2$, the initial value is $u(x, y, 0) = 1$ in $x^2+y^2 < 0.5$, and otherwise $u(x, y, 0) = 0$. The source component is $f = 0$. The homogeneous boundary condition is also assumed. The parameters are selected as $\varepsilon = 0.01$, $T = 0.5$.

This problem also permitted discontinuous initial data and was first given in ^[8] by Karlsen and Risebroand, and since then was widely used for experimenting several numerical methods for nonlinear convection diffusion equations, see ^{[7,9,12,13,14]} and also ^[26,27] for slightly different problem domain or terminal time setting.

For this problem, we know that $u = 0$ for any $(x, y)\in \partial \Omega$ and $t\ge 0$, and by direct calculation

hence, $F_u(u) = (P_u(u), Q_u(u))^T = (0, 0)^T$, and $F_u(u(x, y, t))\cdot{{\bf n}} = 0$, so the inflow boundary $\Gamma_n^{-}$ is always empty for the convection subsystem (2.11) at each time step $t = t^n < T$.

Now we execute Algorithm A under the discretization parameters $\Delta t = 0.001$ and $h = 1/128$ up to the final time $T = 0.5$. Firstly, the contour plot of the numerical solution at $T = 0.5$ is displayed in Figure 5, where the contour values with Labels 1–19 correspond to the lines located from outside to inside in order.

The three dimensional plot of the computational numerical approximations at different times are drawn in Figure 6. From these figures, we can observe the smooth evolution of the numerical solutions by Algorithm A, which verifies the stability of the proposed algorithm. In summary, the effectivity of the present method is verified by good coincide of this contour with those ones in ^{[7,9,12,13,14]}.

4.   Conclusions

We have introduced a new IMEX approach for computing the unsteady nonlinear convection-diffusion problems.

Convection and diffusion subsystems can be both dealt with independently, and moreover, the stiffness matrix of the subsystems keep invariant along time marching. We can use many efficient preconditioned iterative solvers to calculate the diffusion subsystems because of its self-adjointness and coerciveness, like preconditioned Conjugate Gradient solver; at the same time, the mass matrices are permitted in convection subsystems, as well as can be explicitly and fast computed when using spatial linear elements. We also provide a multistep technique to relax the instability, that mostly results from the explicit operation of the convection equations. All the above make the present approach a fast scheme for computing the unsteady nonlinear systems. In the end, we carry out a number of numerical simulations to evaluate the stability and validation of the algorithm. As for the spatial finite element approximation, some special triangulations like the Shishkin-type meshes ^[28] and Bakhvalov meshes ^[29,30], or adaptive grid method ^[31] along time marching can be incorporated with the present time-splitting scheme for the small diffusion coefficient in the future. Also the extensions of the present method to the three-dimensional problems (with complex domains) are possible and will be our further study.

Acknowledgments

The second and third authors are supported by Shenzhen Technology Project grant GXWD20201230155427003-20200822102539001. The fourth author is supported by Fundamental Research Funds for the Central Universities grant 2232019D3-39.

Conflict of interest

The authors declare there is no conflicts of interest.