Numerical modelling of convection-diffusion-reaction problems with free boundary in 1D

TL;DR

This paper introduces a numerical method for solving 1D convection-diffusion-reaction problems with free boundaries, transforming the problem into a fixed domain and modeling interface evolution via ODEs, applicable to various free boundary problems.

Contribution

It presents a novel numerical approach that models interface evolution using local wave shapes and transforms free boundary problems into fixed domain problems for efficient computation.

Findings

The method accurately approximates problems with finite interface speed.

It can handle some problems with infinite interface speed.

Numerical results agree well with analytical solutions.

Abstract

We discuss a numerical method for convection-diffusion-reaction problems with a free boundary in 1D. The method is based on the numerical modelling of the interface evolution, the transformation to a fixed domain problem and the approximation by an ODE system. The interface evolution is modelled by means of the local shape of the corresponding travelling wave solution. The method can be applied to many free boundary problems with a finite speed of the interface. The presented method can also approximate some problems with an infinite speed of the interface for damped travelling wave type solutions. In the numerical experiments we compare our numerical solution with the analytical ones for some problems.

Figures

Tables

Table 1. Table 1: N 𝑁 N versus A ​ L 𝐴 𝐿 AL

$N$	$AL{.10}^4$
10	2.66
20	2.45
30	2.35
40	2.1
50	2.2
100	2.3

Table 2. Table 3: Optimal parameters versus A ​ L 𝐴 𝐿 AL

$N$	$D_p$	$M$	$d$	$AL{.10}^4$
100	$D_3$	68	3	1.2
100	$D_2$	90	3	1.1
100	$D_1$	85	6	1.3
55	$D_3$	16	3	0.55
50	$D_2$	5	5	0.41
50	$D_1$	43	4	1.71

Equations

$${\partial }_{t}u={\partial }_x^2{u}^n+b_0{\partial }_x{u}^{\gamma }+c_0{u}^m$$

$$a)u=\varphi (t),\text{mbox}orb)-{\partial }_x{u}^n+b_0{u}^{\gamma }=q(t).$$

$$u(x,0)=u_0(x)\ge 0,$$

$${\partial }_{t}u={\partial }_x^2{u}^n,n>1,x\in (-\infty ,+\infty )$$

$$u(x,t)=\left\{\begin{array}[a]{c c}\frac{1}{s(t)}[1-(\frac{x}{s(t)})^{2}]^{1/(n-1)},&\mbox{for}\ |x|<s(t)\\ 0,&\mbox{for}\ |x|\geq s(t)\end{array}\right.$$

$$s(t)=[\frac{2n(n+1)}{n-1}(t+1){]}^{1/(n+1)}$$

$$u(x,0)=\left\{\begin{array}[a]{c c}\frac{1}{s(0)}[1-(\frac{x}{s(0)})^{2}]^{1/(n-1)},&\mbox{for}\ |x|<s(0)\\ 0,&\mbox{for}\ |x|\geq s(0).\end{array}\right.$$

$$\dot{s}(t)=-\frac{n}{n-1}{\partial }_x{u}^{n-1}{|}_{x\nearrow s(t)},$$

$${\partial }_{t}u={\partial }_x^{2}a(u)+{\partial }_{x}b(u)+c(u)$$

$${\bar{u}}_t=\bar{\mathcal{L}}(\bar{u})-\frac{\dot{s}}{s}y{\partial }_y\bar{u},y\in (0,1)$$

$$\dot{s}=\mathcal{B}({\partial }_y\bar{u}(1,t)),\text{mbox}resp.\dot{s}=\mathcal{B}({\partial }_y\bar{u}(1,t),{\partial }_y^2\bar{u}(1,t)),$$

$$\left\{\begin{array}[a]{c c}f(\xi)=0,&\mbox{for}\ \xi\geq\xi_{0}\\ f(\xi)>0,&\mbox{for}\ \xi<\xi_{0}\end{array}\right.$$

$$-\sigma f^{\prime }=n{f}^{n-1}{f}^{\prime \prime }+n(n-1)f^{n-2}(f^{\prime }{)}^2+b_0\gamma f^{\gamma -1}(f){)}^{\prime }+c_{0}f,\xi \in (-\infty ,\infty ).$$

$$\theta (f(\xi ))=-(f^n{)}^{\prime }(\xi ),$$

$$\theta (s)=\sigma s+b_0{s}^{\gamma }-n{c}_0{\int }_0^s\frac{r^{m+n-1}}{\theta (r)}dr.$$

$$\theta (s)\approx {\theta }_0{s}^q\text{mbox}ass\searrow 0$$

$$\begin{array}[a]{c c c c}(i)&c_{0}<0;&\quad(\mbox{in this case}&q=\min\{(n+m)/2,1\})\\ (ii)&c_{0}=0;&\quad(\mbox{in this case}&q=1)\\ (iii)&c_{0}>0,\ m+n\geq 2;&\quad(\mbox{in this case}&q=1).\end{array}$$

$$\theta (s)={\theta }_0{s}^q+o(s^q),\text{mbox}with{s}^{-q}o(s^q)\to 0,\text{mbox}fors\searrow 0.$$

$$f(\xi )\approx d^{\beta }({\xi }_0-\xi {)}^{\beta }\text{mbox}for\xi \nearrow {\xi }_0.$$

$$f^{\prime }(\xi )\approx \frac{{\theta }_0}{n}f(\xi {)}^{\alpha },\alpha =q-n+1$$

$$f(\xi )\approx [\frac{{\theta }_0(n-q)}{n}{]}^{\beta }{(\xi }_0-\xi {)}^{\beta },\text{mbox}where\beta =\frac{1}{n-q},d=\frac{{\theta }_0(n-q)}{n}$$

$$f(\xi )\approx d^{\beta }({\xi }_0-\xi {)}^{\beta },f^{\prime }(\xi )\approx \beta d^{\beta }({\xi }_0-\xi {)}^{\beta -1}\text{mbox}for\xi \nearrow {\xi }_0,$$

$$[f^n{]}^{\prime }(\xi )\approx d^{n\beta }n\beta ({\xi }_0-\xi {)}^{n\beta -1}\text{mbox}for\xi \nearrow {\xi }_0,$$

$$[f^n{]}^{\prime \prime }(\xi )\approx d^{n\beta }n\beta (n\beta -1)({\xi }_0-\xi {)}^{n\beta -2}\text{mbox}for\xi \nearrow {\xi }_0$$

$$\dot{s}(t)=Q(w)=\left\{\begin{array}[a]{c c}-\frac{n}{(n-1)}\partial_{x}w+\frac{c_{0}}{\partial_{x}w}-G(\gamma)b_{0}&\ \mbox{for}\ m+n=2;\\ -\frac{n}{(n-1)}\partial_{x}w-G(\gamma)b_{0}&\ \mbox{for}\ m+n>2\\ -\frac{n}{(n-1)}\partial_{x}w-G(\gamma)b_{0}&\ \mbox{if}\ c_{0}=0\par\end{array}\ ,\quad x=s(t)\right.$$

$$\dot{s}(t)=-\frac{{\partial }_{t}u}{{\partial }_{x}u}\text{mbox}inthesensex\nearrow s(t),$$

$$\dot{s}(t)=\underset{x\nearrow s}{\lim }\frac{1}{d^{\beta }(s-x{)}^{\beta -1}+o_1}\{d^{n\beta }n\beta (n\beta -1)(s(t)-x{)}^{n\beta -2}-$$

$$b_0{d}^{\gamma \beta }\gamma \beta (s-x{)}^{\gamma \beta -1}+c_0{d}^{m\beta }(s-x{)}^{m\beta }+o_2+o_3+o_4\}$$

$$o_0\equiv o((s-x{)}^{\beta }),o_1\equiv o((s-x{)}^{\beta -1}),o_2\equiv o((s-x{)}^{n\beta -2}),$$

$$o_3\equiv o((s-x{)}^{\gamma \beta -1}),o_4\equiv o((s-x{)}^{m\beta }).$$

Full text

Keywords: interface modelling, free boundary problems, porous media equations, degenerate diffusion-adsorption

1 Introduction

We shall consider convection-diffusion-reaction equations of the type

[TABLE]

on the halfline $x\in(0,\infty)$ for $t>0$, with the boundary conditions

[TABLE]

at the boundary $x=0$. Here, $b_{0}$ and $c_{0}$ are constants. The initial condition is of the form

[TABLE]

where $u_{0}(x)>0$ for $x\in(0,L_{0})$ and $u_{0}(x)=0$ for $x>L_{0}$.  We are only looking for nonnegative solutions since negative solutions have no physical interpretation. We assume that the parameters $n,\gamma$ and $m$ are nonnegative real numbers and guarantee the existence of a function $s(t)$ such that the solution is positive in $(0,s(t))$ and $u(x,t)=0$ for $x>s(t)$. The function $s(t)$ represents the time evolution of the interface. Such solutions have been intensively studied in the last 30-years in many papers, e.g., O.A. Olejnik, A.S. Kalashnikov, Czou-Ju-Lin [22], G.J.Barenblatt [3], D.G. Aronson [1], B.H.Gilding [11, 12], L.A.Peletier [23], R.Kersner [21], B.H.Gilding and R.Kersner [13, 14] etc. (see also the list of references in the previous papers).

We shall discuss the numerical approximation for (1)-(3) based on interface modelling. In these problems the free boundary $s(t)$ is not given explicitly in the model and its determination is included in the solution. Solutions to the linear convection-diffusion-adsorption problems ($n=m=\gamma=1$) in (1) do not show the interface property. There, an initially localized support of the initial profile spreads with an infinite speed ($\dot{s}(t)=\infty$). Degeneration of the elliptic part ($n>1$) is the main reason for the appearance of the interface. However, also in the regular case ($n=1$) the interface will appear under special properties of the adsorption represented by $c_{0}$ and $m$. Very substantial results (necessary and sufficient conditions for $n,m,\gamma$ in (1)) to guarantee the appearance of the interface have been obtained in many papers, e.g. B.Song ([25]), B.H.Gilding and R.Kersner ([13, 14]), see also list of references there. If $b_{0}\neq 0,c_{0}<0$ and $m,n,\gamma>0$ the existence of an interface is guaranteed when $\min(n,\gamma)>\min(m,1)$. This condition is also necessary (see [25]). If $b_{0}=0$ and $c_{0}=0$ in (1) then the condition $n>1$ is the real criterion whether or not the interface appears. An equation with this property is called “porous media equation” and an analytical solution has been found by Barenblatt [3]. This equation is a mathematical model for adiabatic infiltration of a gas into a porous media. Its solution is not a classical one (but a generalized, weak solution, defined by a corresponding integral identity) since $\partial_{x}u$ is not continuous along the interface, i.e. $x=s(t)$. The Barenblatt-Pattle solution of the Cauchy problem for the porous media equation

[TABLE]

is of the form

[TABLE]

where

[TABLE]

and the corresponding initial condition is

[TABLE]

This analytical solution stimulated the theory of existence and uniqueness for the more general “porous media” type problems. More detailed information on this topic can be found in [13, 14] and the references therein.

The solution of (4) for $n=6$ and $t\in(0,200)$ is drawn in Fig.1 (in $10$ equidistant time sections) and the evolution of its interface is drawn in Fig.2. We can readily verify the validity of the equation

[TABLE]

which is the governing ODE model for this iterface. Solutions of the problem (1)-(3) have the shape of damped travelling waves with eventually sharp fronts at the interface. This makes the numerical approximation very difficult (the movement of the interface is not priori known). There exist numerical methods for free boundary problems that estimat the position of the interface and then apply a local adaptation of the space discretization, which are very expensive but have a limited precision (see e.g. [17, 18] etc.). Some other results related to our method can be found in [4, 6].

The main goal of this paper is the construction of an efficient and precise numerical method for solving “porous media” type problems (1)-(3) . This method is based on the modelling (in terms of an ODE) of the time evolution of the interface $s(t)$ , which significantly depends on the local profile of the solution near the interface. The unknown $s$ will be explicitly included into the solution of the problem. Then, we shall look for the solution in the interval $(0,s(t))$ only (instead of $(0,L)$, where $L>L_{0}$). Consequently, we will transform (1) to the fixed domain $(0,1)$ (using $y=\frac{x}{s(t)}$ ) and invoke a special space discretization (finer in the neighbourhood of the point $y=1$ which corresponds to the interface $x=s(t)$). This approximation of (1) results in a corresponding system of ODEs including the uknown $\dot{s}(t)$. This system is completed with the ODE model for the interface (of the same type as (6) for (4)). The final ODE system is solved with an appropriate solver for stiff ODE systems. Thus, the modelling of the interface evolution is the crucial point for the construction of our numerical method. In this way we can approximate the sharp front of the solution at the interface with higher accuracy. In our numerical experiments we have used the solver ’ode15s’ from MATLAB library.

In Section 1 we describe the method in more detail. In Section 2 we construct the interface model for (1) under some assumptions on parameters $n,m$ and $\gamma$. In section 4 we construct the numerical approximation for (1). In Section 5 we extend the results from the previous sections to the more general PDE

[TABLE]

satisfying $a(0)=b(0)=c(0)=0$, $a^{\prime}(0)=0$ and preserving similar interface property of the solution as (1) . Also more general mathematical models involving $a(t,u),b(t,x,u),c(t,x,u)$ can be included too. These types of equations describe a large part of convection-diffusion-adsorption problems with many important practical applications. In Section 6 we discuss the interface modelling in the non degenerate case $a^{\prime}(0)>0$. These equations are of the “oxygen-consumption” type, where appearence of the interface is generated by strong adsorption represented by $c(u)$. In Section 7 we present some numerical experiments and we present comparisons with the analytical solutions.

2 Idea of the numerical approximation

For a good numerical approximation of (1) it would be desirable to use space discretization only in the moving domain $(0,s(t))$ and to use also moving grid points which follow the interface and which are more dense at the sharp fronts near the interface - see Fig.1. For this purpose we formally consider (1) in the domain $(0,s(t))$ which by Landau’s transformation, $y=\frac{x}{s(t)}$, is transformed it to the fixed domain $(0,1)$. Then, for $\bar{u}(y,t):=u(x,t)$ we obtain the equation

[TABLE]

where ${\cal L}(u)$ is the rhs of (1) and $\bar{{\cal L}}(\bar{u})$ is its transform in the $y$-variable. If we can model $\dot{s}$ in terms of $\partial_{y}\bar{u}(1,t)$, resp. $\partial_{y}^{2}\bar{u}(1,t)$ (which describe local properties of $u$ in a ’small’ neighbourhood of $s(t)$), e.g., in the form (see (1))

[TABLE]

then we can eliminate $\dot{s}$ in (8) by means of the rhs of (9). Then the complete system consists from (8) and (9). We approximate it using the space discretization $\{y_{i}\}_{i=1}^{N}\in(0,1)$ of the $y$-variable. The resulting system of ODEs is stiff and we solve it by a corresponding ODE solver. Now, the unknown interface is explicitly included in the solution. The grid points $\{y_{i}\}_{i=1}^{N}$ correspond to the moving grid points $\{x_{i}(t)\}_{i=1}^{N}\in(0,s(t))$, with $x_{N}(t)=s(t)$ for (1).

The method used cannot be extended to more space dimensions (except of radial symmetric ones). The mathematical model for the interface movement (its normal velocity) depends not only on local properties of the solution (e.g., the normal space derivative) but also on the geometry of the support of the solution, which in 1D is very simple.

3 Modelling of the interface

The appearance of the interface in (1) is closely related to the existence of the travelling wave solution (see [13, 14] ) of the form $u(x,t)=f(\xi),\ \xi=x-\sigma t$ where

[TABLE]

for some $\xi_{0}\in(-\infty,+\infty)$. The function $f$ has to be determined from the ODE

[TABLE]

The function $\theta$ defined as

[TABLE]

has to satisfy the nonlinear Volterra equation (see [13])

[TABLE]

In [13] (Lemma 13) one proves: There exists a solution $\theta$ of (10) for all $\sigma>\sigma^{\ast}$ (for suitable $\sigma^{\ast}>0$) provided the conditions (i)-(iii) below with parameters $n\geq 1,\ m>-n$ and $\gamma\geq 1$ are fulfilled. Moreover, there exists a $\sigma^{\ast\ast}>\sigma^{\ast}$ so that

[TABLE]

for some $\theta_{0}$ and $q>0$ depending on $n,m,\gamma$. These conditions are

[TABLE]

The conditions (i)-(iii) are also necessary for the existence of the solution of (10).

The case $0<n<1$ is also discussed in [13] (Lemma 13), but we cannot aply it in our analysis.

The relation (11) is equivalent to

[TABLE]

As a consequence of (11), also the front of the semi-wave $f$ (localized at the point $\xi_{0}$), satisfies the approximation (see Remark 1 below)

[TABLE]

We shall use this property in the modelling of $\dot{s}(t)$ in terms of local properties of the solution at the interface (e.g. $\partial_{x}u(s(t),t)$) and model parameters. This will be substantially used in our numerical method. In fact, we shall assume the asymptotic property (12) to hold also for the first and second derivatives of $f$ (see (13) below).

Remark 1. From (11)and $\theta(f(\xi)):=\frac{d}{d\xi}f^{n}(\xi))$ (provided $f$ is smooth, decreasing and $f(\xi)=0$ for $\xi\geq{\xi}_{0}$ - “travelling semi-wave”) it follows

[TABLE]

and consequently, by integration over $(\xi,{\xi}_{0})$, we obtain (see (12)),

[TABLE]

under the assumption $n\neq q$.

The development of an interface model of the type (9) in a rigorous mathematical way requires a very deep analysis. Moreover, there are no results available in this respect (except of a few, very special cases). We will construct such a model under very strong smoothness assumptions on the solution in the interval $(0,s(t)]$ (up to the free boundary) and the assumption (12) concerning the shape of travelling semiwave corresponding to (1). In fact we assume (12) in its stronger form

[TABLE]

[TABLE]

[TABLE]

where $\beta$ and $\ d$ are unknown parameters to be determined.

Our main result in this section is

Theorem 1. Suppose (13) and $n>1,\ 0\leq m<1,\gamma\geq 1$ . When the solution of (1) is regular up to the interface (the equation is satisfied in the limit sense at the interface), then the interface $s$ of (1) is governed by ODE

[TABLE]

where $\beta=\frac{1}{n-1}$,  $u=w^{\beta}$ and $G(r)=0$ for $r>1$, $G(1)=1$.

Proof. Our starting point is the condition $u(s(t),t)=0$ for all $t\in(0,T)$. Differentiating this equation with respect to $t$ we obtain formally

[TABLE]

which cannot be used practically, since usually we have $\partial_{x}u\to-\infty$ for $x\nearrow s(t)$. But invoking the assumption (13), we can conclude that $\partial_{x}u^{\frac{1}{\beta}}$ is finite for $x\nearrow s(t)$ for some $\beta>0$ . Using the smoothness of the solution (the equation (1) holds up to the interface $x=s(t)$) we can eliminate $\partial_{t}u$ in (15) by the rhs of (1) and then we insert $u=d^{\beta}(s(t)-x)^{\beta}+o_{0}$ in the small neighbourhood of the interface (for $x\nearrow s(t)$). Due to (13), we obtain

[TABLE]

[TABLE]

with the lower order terms represented by $o_{1}-o_{4}$ ( because of (13)), where

[TABLE]

[TABLE]

The assumtions on $n,m$ and $\gamma$ guarantee the existence of the interface, i.e. $|\dot{s}(t)|<\infty$ because of [13] (Lemma 13). Our goal is to find such a $\beta>0$ that this will be achieved. There are more combinations of parameters $n,m,\gamma$ to guarantee that such a $\beta$ could be found. Deviding the nominator on the rhs in (16) by $(s-x)^{\beta-1}$, we obtain three terms containing $s-x$ with the exponents $(n-1)\beta-1;\ (\gamma-1)\beta;\ (m-1)\beta+1$ (corresponing to diffusion, convection and reactive part) . To guarantee $|\dot{s}(t)|<\infty$, we have to choose $\beta>0$ such that all three exponents are nonnegative and at least one of them has to be zero. The case $\beta\leq 0$ doesn’t correspond to the required property of the semiwave. Then, for $n+m\geq 2$ the choice $\beta=\frac{1}{n-1}$ fulfilles our requirements. This choice also includes the case $c_{0}=0$ . Then from (16) we deduce the formula for $\dot{s}(t)$ which contains the uknown parameter $d$ -see (12). In the following we shall determine $d$ . For this purpose consider the new unknown $w$ defined by $u=w^{\beta}$ and rewrite (1) in terms of $w$. In the small neighbourhood of the interface ($x\nearrow s(t)$) we have $w(x,t)\approx d(s(t)-x)$ due to (13) and hence

[TABLE]

From this we conclude that $\partial_{x}w\to-d$ when $x\nearrow s(t)$.

Now, we substitute $d=-\partial_{x}w$ into (16) and we obtain the required ODE model for the interface.

Remark 2. It is very difficult to verify the assumptions (13). However consequences from them, leading to the interface models, are verified in some examples in Section 7, where the analytical solution is available. Since assumptions (13) could be violated in some model problems, our modelling of the interface should be justified numerically in practical applications.

Remark 3. The analysis used in choosing $\beta$ so that $|\dot{s}(t)|<\infty$ is rough and cannot distinguish the sign of $b_{0}$ and $c_{0}$ which is important for the existence of the interface. The analysis for the adsorption and reaction is the same. From (i) and (iii) it follows that the interface exists also in the case $n+m<2$ provided $c_{0}<0$ (i.e. in the adsorption case), while for $c_{0}>0$ the interface doesn’t exist.

Remark 4. In the case $n+m\geq 2$ it follows from (i),(iii) that $q=1$. Then, our $\beta$ obtained in the proof of Theorem 1 coincides with the one from Remark 1 what supports our arguments. In fact, the assumptions in Theorem 1 include also the additional solution $\beta=\frac{1}{1-m}$ in the case $n+m>2$. However, this $\beta$ is excluded by the fact $q=1$ (see (i), (iii) and Remark 1). It means that our analysis leading to Theorem 1 is not satisfactory for the unique determination of $\beta$.

4 Numerical approximation of (1)

The mathematical model for the time evolution of $s$ given by Theorem 1 is of practical use, since $\partial_{x}w$ is finite and nonzero. From this reason we also rewrite (1) in terms of $w$ using the transformation $u=w^{\beta}$ where $\beta=\frac{1}{n-1}$ . Then, we obtain

[TABLE]

We have to solve this equation in the moving domain $(0,s(t))$. We transform (17) to the fixed domain $(0,1)$ using Landau’s transformation $y=\frac{x}{s(t)}$, with $\bar{w}(y,t)=w(x,t)$. For the sake of simplicity, in the following we write $w$ in the place of $\bar{w}$. Then, we obtain

[TABLE]

[TABLE]

This equation has to be completed with the ODE for $\dot{s}$ given by Theorem 1 and $\dot{s}$ at the rhs has to be eliminated by means of the rhs of (14) . Then, we have to solve the coupled PDE and ODE system

[TABLE]

[TABLE]

[TABLE]

on a fixed domain $y\in(0,1),\ t>0$, where $Q$ is from (14) (there, $\partial_{x}w$ has to be replaced by $\frac{\partial_{y}w}{s}$). The PDE (4) at the point $(y,t)$ depends also on the point $(1,t)$ because of $Q(w)$. Consequently, it is non-local. In the numerical approximation of (4)-(19) we use space discretization (method of lines) which results in an ODE system. We generally use a nonequidistant partition of $(0,1)$. Let $\alpha_{i}>0,\ \forall i=1,...,N$,  $\alpha_{0}=0$ and consider the grid points $y_{i}=\sum_{j=0}^{i}\alpha_{j},\ \forall i=0,...,N$. We approximate $w(y_{i},t)\approx C_{i}(t)$. Denote by $\ell_{i}(y)$ the Lagrange polynomial of the second order through the points $(y_{i-1},C_{i-1})$, $(y_{i},C_{i})$ and $(y_{i+1},C_{i+1})$. We approximate $\partial_{y}^{2}w(y_{i},t)$ and $\partial_{y}w(y_{i},t)$ by the corresponding expressions given by $\frac{d^{2}\ell_{i}}{dy^{2}}|_{y=y_{i}}$ and $\frac{d\ell_{i}}{dy}|_{y=y_{i}}$.

We approximate the derivative $\partial_{y}w(1,t)$ appearing in $Q$ by

[TABLE]

where $\ell_{B}$ is the second order Lagrange polynomial throudh the points $(y_{N-2},C_{N-2})$ and $(y_{N-1},C_{N-1})$, $(1,0)$. Our ODE system now reads

[TABLE]

[TABLE]

[TABLE]

for $i=1,...,N-1$, with $\beta=\frac{1}{n-1}$. In the case of Dirichlet boundary conditions this is the resulting system of ODE, where we use $C_{0}(t)\equiv\phi(t)$ (see (2), case a)). In the case of a Robin condition (see (2), case b)), the system (20)-(21) has to by completed by the additional ODE at the point $y=0$. In this case we make use of a ghost point $y_{-1}=-\alpha_{1}$. The missing value $C_{-1}$ is determined from the boundary condition $(\ref{2})$, (case b), by means of $q$ and $C_{1}$ from the relation

[TABLE]

Then, the additional ODE in $y=0$ is the same as (20) for $i=0$ using $C_{-1}$.

In some physical problems it is important to know the support of the solution which in 1D is given by the position of the interface. This is explicitly included in our approximation represented by (20)- (21). This gives us the possibility to determine $s(t)$ with higher accuracy then the approximations without interface modelling.

Remark 5. The system (20)-(21) is stiff and a corresponding ODE solver has to be used. This system is a good starting point to solve the original free boundary value problem. The first advantage is that we approximate the solution on its support only. The second one lies in the space discretization. In fact, we have moving grid points (in the original $x$- variable) by means of which sharp fronts near the interface can be approximated very well for all $t$. We shall demonstrate this in our numerical experiments presented in Section 7 . The presented numerical approximation requires positive initial profile on its support $(0,s(0))$. If our initial profile is zero (provided nonzero boundary condition at $x=0$ is considered), then we have to start with the “small” positive auxiliary initial profile on the “ small ” domain $(0,s(0))$. Also some compatibility between the initial and boundary conditions (at the point $(t,x)=(0,s(0))$) help to avoid numerical instabilities in the ODE solver at the starting time point. This “small” regularization does not change the expected solution significantly. The efficiency of the numerical method considerably depends on the strategy of grid points distribution. This allows reducing the size of the system (20) without decreasing the accuracy.

The used space discretization of the PDF and its approximation it by an ODE- system is of the second order and the order of ODE solver is regulated by its tolerance parameters.

5 Extension of the method to PDE (7)

From the theory of porous media type solutions one can notice that the evolution of the front of the solution substantially depends on the order of degeneration and adsorption and on local properties of the solution. Therefore, the same idea of the numerical approximation used on (1) can also be applied in a more general case of (7), where we replace $a(s)\leftrightarrow a(t,s),\ b(s)\leftrightarrow b(x,t,s),\ c(s)\leftrightarrow c(x,t,s)$, provided $a,b$ and $c$ are asymptotically polynomial with respect to $s$ at the point $s=0$. More precisely, we assume

[TABLE]

[TABLE]

Modelling the interface in this setting (generalization of (7)), we proceed similarly as in Section 2 and obtain.

Theorem 2. Let the assumptions of Theorem 1 hold true. Moreover, suppose (22). Then, the governing ODE for the interface in the genaralized case of (7) is given by

[TABLE]

where $u=w^{\beta}$ with $\beta=\frac{1}{n-1}$ and $G(r)=0$ for $r>1$, $G(1)=1$ and $u=w^{\beta}$ with $\beta=\frac{1}{n-1}$.

Using this ODE model for the interface, we approximate the considered PDE similarly as in the case of (1).

6 Modelling of the interface in the nondegenerated case $\partial_{u}a(t,0)>0$.

Cosider the generalization of (7) as in the previous section with the nondegenerate elliptic part $\partial_{u}a(t,0)>0$. As it was mentioned in the introduction, the interface can appear also in the nondegenerated case when adsorption is sufficiently strong for small values of $u$ . As an example we can take oxygen-consumption type problem

[TABLE]

with $u(x,0)=u_{0}(x)>0$ for $x\in(0,L_{0})$, $u_{0}(x)=0$ for $|x|>L_{0}$, where $D$ is the diffusion coefficient and

[TABLE]

represents the consumption of the oxygen with concentration $u$. This model problem has also an interface with finite speed and has been analysed some decades ago (see, e.g., [9]). There exists an analytical solution wich is a polynomial of the second order at the front (near the interface). This motivates us to look for a polynomial solution at the interface also in the more general case (7) considered in the previous section. The order of this polynomial can be determined from the information that the considered problem has a finite speed of interface propagation. Then, for the interface modelling we proceed similarly as in the degenerated case. The existence of the interface (in case of (7)) has been established in [13] (Theorem 12, Assertion (i)).

For simplicity we consider a generalization of the mathematical model (7) (of oxygen-diffusion type) under the assumptions

[TABLE]

[TABLE]

A more special case when $c=c_{0}H(u)$ has been treated numerically in [7] and we have extend this idea to the more general case of (7), where (25) is fulfilled. Another approximation method based on the solution of variational inequalities has been used in [3].

We start with the formal “interface condition” $\dot{s}=-\frac{\partial_{t}u}{\partial_{x}u},$ where we substitute the rhs from (7) with $\partial_{t}u$ under the smoothness assumption on the solution and with the validity of (7) up to the interface. Then, we get

[TABLE]

where the rhs is understood in the limit sense $x\nearrow s(t)$. Again, we consider the transformation $u=w^{\alpha}$ for some $\alpha>0$ in (7) in the setting (25). We obtain

[TABLE]

where the variables $x$ and $t$ in functions $a_{0},b_{0}$ and $c_{0}$ are omitted. We are looking for such an $\alpha$ that (27) holds up to the interface $x=s(t)$. Since $w(s(t),t)=0$ ($w(x,t)\to 0$ for $x\to s(t)$), it must hold that

[TABLE]

Due to $a_{0}(s(t),t)\geq\delta>0$, we conclude

[TABLE]

The negative sign is chosen due to the fact that $w$ is decreasing as $x\nearrow s(t)$. Since

[TABLE]

and $|\dot{s}(t)|<\infty$, we choose $\alpha$ so that $\partial_{x}w\neq 0$. Then, with respect to (29) we choose

[TABLE]

This guarantees that $|\partial_{x}w|<\infty$. To make use of (27) in (30) we have to compute the limit

[TABLE]

Both, the nominator and denominator tend to zero as $x\nearrow s(t)$, so the L’Hospital rule can be used. We obtain

[TABLE]

[TABLE]

[TABLE]

Finally, from (26)-(31) we conclude

[TABLE]

[TABLE]

Theorem 3. Consider (7) under the assumption (25). Let the solution $u$ be $C^{2}$- smooth up to the interface $s(t)$. Then, the solution admits an interface $s(t)$ and this is governed by the ODE (32).

In the case $m=0$, the solution $u$ is asymptotically a quadratic polynomial at the interface, since $u=w^{2}$, $\partial_{x}w\neq 0$, and $u$ is finite on the interface. This is the case for a classical oxygen-consumption problem (24).

This result can be extended to the more general case (7), with $a=a_{0}(x,t)g(u)u,b=b_{0}u^{\gamma}$, where $g(0)>\delta>0,\ \gamma\geq 1$ and $c$ is as in (25).

The numerical approximation of the generalized oxygen-consumption problem (7) under the conditions (25) proceeds along the same lines as in Section 4, using our governing equation (32). The approximation $\partial_{x}^{2}w$ on the interface by means of $\partial_{y}^{2}\ell_{B}|_{y=1}$ is generaly not satisfactory. The third order polynomial through the points

[TABLE]

has to be used.

7 Numerical experiments

To make an advantage of our numerical approximation a suitable discretization (in the $y\in(0,1)$-variable) can be chosen. Since $y=1$ corresponds to the interface, where the front of the travelling and damping semiwave is situated, we increase the density of grid points in the neighbourhood of $y=1$. At the same time we decrease the density of grid points at the trailing part of the damping semiwave type solution to keep the set of all grid points small. This of course depends on the solution profile and on the type of the solved model. The distribution of the grid points plays an important role in the accuracy and efficiency of the numerical realization. Numerical approximation without the interface information requires much more grid points in the space discretization then our method with interface modelling to obtain the same order of accuracy. For large $t$ the profile of the solution can be significantly changed comparing with that one for small $t$. In that case we can restart the numerical solution of the corresponding ODE-system for large $t$ with a new grid (with a changed strategy of grid points distribution) and the corresponding initial profile. The main goal is to obtain the same required accuracy with a minimal set of grid points, which determines the dimension of our ODE system to be solved. We have used the following discretizations. The interval $(0,1)$ is divided equidistantly in $M$ subintervals and the last $d$ intervals are successively subdivided :to $2^{p},3^{p},...(d-1)^{p}$ additional subintervals for $p=1,2,3$. We denote these discretizations by $D_{p}$. Mostly, we will apply a smooth geometrical dicretization, where the length of the $(i+1)$-th subinterval is a $q<1$ factor smaller than that of the $i$-th. At this strategy, $q$ is determined from the prescribed number of all grid points $N$ and the length of the first subinterval $1/M$. We denote this discretization by $D_{4}$. Space discretization $D_{4}$ seems to be most efficient in most of our numerical applications.

In our numerical experiments we solve the resulting ODE-system by means of solver ode15s from MATLAB library, which is developed for stiff ODE-systems.

7.1 Barenblatt solution.

The interface model given by Theorem 1 (the case with $c_{0}=0$) coincides with that one in (6). In Fig. 1 we present our numerical solution using $n=6$ and the discretization $D_{4}$ with parameters $M=10$ and $N=20$. The corresponding analytical solution cannot be distinguished graphically from the numerical one. The time evolution of the interfaces is shown in Fig. 2. The interfaces cannot be distinguished, too, even during the long time period $(0,200)$.

To compare the numerical $u_{n}$ and analytical $u_{an}$ solutions, we use a (relative) error estimator given by a numerical $L_{2}$-norm . We use discrete time points $t_{j}$ of the solutions approximated in $x_{i}=s(t_{j})y_{i},\ i=1,...,N$ grid points and define

[TABLE]

where the space integration is performed over the maximum length of supports of both solutions. In the same way we define the $L_{1,rel}(t_{j})$ error (the exponent $2$ is replaced by $1$ in the integrands). Mostly (in our experiments below) the time evolution of both errors does not differs significantly. The time evolution of the relative $L_{2}$- error is shown in Fig. 3.

Using more grid points ($N=40,M=15$), the maximum $L_{2}$-error is $2.2e-04$ and for ($N=60,M=20$) it is $1.3e-04$.

Let us approximate (4) in a similar way (reducing it to an ODE-system) without the interface modelling. We use a mass conservation scheme where

[TABLE]

with $\alpha(i)=y_{i}-y_{i-1},\ i=1,...,N-1$. At the point $y=0$ we obtain the ODE

[TABLE]

due to the symmetry arguments ($\partial_{x}u|_{x=0}=0$). In Fig. 4 we present the corresponding numerical solution for $p=6$ with an equidistant space discretization of $x\in[0,10]$ with $N=M=100$ grid points. The error $L_{2,rel}(t)\in(0.13,0.018)$ for $t\in(0,200)$. This error is nearly 100 times higher than the error of numerical results based on interface modelling with only $N=20$ grid points. The relative $L_{2,rel}(t)$-error for the classical numerical approximation increases significantly with the reduction of grid points. The most of the contribution to the error is located near the interface.

In the following Tables 1-3 we present the relation between the space discretization $D_{p},\ p=1,2,3,4$ and the average of $L_{2,rel}(t)$-error versus the time $t\in(0,200)$ in terms of the following numerical approximation

[TABLE]

where we take $r=30$. We compare our numerical solution (based on interface modelling) with the analytical Barenblatt-Pattle solution for $n=6$. In Table 2 we present $N$ and $AL$ where $M=N$ and discretization $D_{4}$ (i.e., unifom discretization in $y\in(0,1)$). In Table 2 we chose an optimal $2\leq M\leq N$ so that $AL$ is minimal. In Table 4 we present $N,D_{p},d,M$ and $AL$ where for given $N$ and $D_{p}$ an optimal $M$ and $d$ are chosen so that $AL$ is minimal.

7.2 An example of porous media equation with adsorption

The time evolution of the interface can be very interesting when degenerated diffusion and adsorption is considered. In the following example we may see the dependence of the time evolution of the interface on the local profile of the solution (at the interface) and observe the order of the degeneracy very transparently. Let us consider the model problem

[TABLE]

with $\partial_{x}u_{x}|_{x=0}=0$. For a special initial profile there exists an analytical solution [21] given by the formula

[TABLE]

[TABLE]

(when the value in parenthesis [.] is nonnegative, otherwise $u=0$) which corresponds to the initial profile

[TABLE]

where $\alpha>0,L_{0}>0$ are given parameters.

This example suits to our model considered in Theorem 1 (where $n=p,m=2-p,\forall p\in(1,2)$). The interface is modelled (see Theorem 1) by

[TABLE]

Since we know an analytical solution (for the special initial profile above) a very important question arises. Does the analytical solution satisfy our interface model given by Theorem 1? The answer is positive and this can be verified by some additional computations. This supports our results in interface modelling. In Fig. 5-8 we present the analytical solution, numerical solution, interfaces and relative error, respectively, for the parameters $p=1.8$, $N=100$ and $M=80$. In this case the interface starts to move forwards in the beginning (there is a sufficiently high slope of the initial profile at the point $L_{0}=s(0)$) and turns backwards at the time $t\approx 5.5$, since then the influence of the adsorption prevails.

At time $t_{que}\in(17,17.5)$ the solution is zero and this is a singularity point for our method. Also in the neighbourhood of $t=17$ the error significantly increases (the solution is very small). In Fig. 9 we present the interfaces using only $N=20$ grid points. In this case $L_{2,rel}(t)\in(0,0.05)$ for $t\in(0,12.5)$ and then increases from $0.05$ up to $0.15$ for $t\in(12.5,17)$.

As we can see even $N=20$ grid points are sufficient for the numerical approximation.

Similarly as in the previous example we solve this problem numerically without the interface modelling. In Fig. 10 we show the time evolution of the numerical solution when using $N=100$ grid points. We consider the domain of the solution to be $(0,20)$.

In this case $L_{2,rel}(t)\in(0,0.05)$ for $t\in(0,12)$ and then increasis from $0.05$ up to $0.25$ for $t\in(12,17)$. When using only $N=40$ grid points, the time evolution of the error is $L_{2,rel}(t)\in(0,0.5)$ for $t\in(0,8)$ and further increases from $0.5$ up to $5$ for $t\in(8,14)$. The front of the solution is not as sharp as in the case of the Barenblat solution. Here, the best discretization strategy is the uniform partition ($N=M$ of the interval $y\in(0,1)$) . The dependence of $AL$ on $N$ is presented in Table 4, where the time interval $(0,14)$ has been considered.

7.3 Convection-diffusion-adsorption model (Contaminant transport)

The governing equation of the (simplified 1D) contaminant transport model with adsorption is of the form

[TABLE]

coupled with the kinetics of adsorption

[TABLE]

where $u$ is the concentration of the contaminant, $v$ is the velocity of the fluid in porous media and $S$ is the adsorbed contaminant by the unit mass of the porous media with the density $\rho$. The function $\Psi$ is an adsorption isotherm, e.g., $\Psi(s)=c_{0}s^{p},p\in(0,1)$ (Freundlich isotherm). When $\kappa\to\infty$ (adsorption is realized in an equilibrium mode), then $S=\Psi(u)$ and in this case (35) is rewritten in the form

[TABLE]

which is of the form (7) after the transformation $\phi=B(u)$.

Numerical modelling of the interface (in the case of Freundlich isotherm) has been discussed in Section 5. For the construction of the interface model we do not invert $B$ and use (37)

[TABLE]

in the interface condition

[TABLE]

In this case the interface is modelled by the ODE

[TABLE]

Implementing this interface model and rewritting (37) in terms of $\phi$, we proceed as in Section 4 and obtain the following numerical results. We use the parameters $p=1/2,\ t\in(0,1.9),\ N=150,M=40,\ \ v=1,\ D=0.05,\ a=\rho=1$ and the Dirichlet condition $u(0,t)=1$ on the boundary $x=0$ . The solution (concentration $u$) is presented in Fig. 11 in 11 equidistant time sections (starting from $t=0$), where the initial concentration is represented by the first (blue) graph (a regularization of the initial zero profile). The corresponding interface evolution is presented in Fig. 12. This will be denoted as the case I. Then, the final concentration profile $u(x,1.9)$ is taken as an initial profile for thecontaminant transport with the homogeneous Dirichlet boundary condition $u(0,t)=0,\ \forall t\in(0,40)$ . The remaining parameters of the model are the same as in the previous experiment. This will be denoted as the case II. The evolution of the concentration profile is presented in Fig. 13 and the corresponding interface in Fig. 14.

The numerical difficulty increases when we consider higher order degeneracy represented by $p=0.25$ and a small diffusion coefficient $D=0.005$. We keep other parameters from the previous experiment. In Fig. 15-16 we present the time evolution ( in $t\in(0,1.9)$-case I and $t\in(0,40)$-case II) of the concentration.

The concentration profile drastically changes in Fig. 16. Therefore, we present the corresponding solution in $11$ time sections of the interval $t\in(0,10)$ (in case II) - see Fig. 17. In this experiment the considered model is very near to nonlinear transport and the solution is very near to the entropy solution of the nonlinear transport with $D=0$

• see [10]. To increase the density of the grid points at the front we take $M=20,N=150$.

In the fact the result presented in Fig. 17 demonstrates the efficiency of our numerical method. The almost piecewise constant initial profile undergoes the nonlinear transport with zero boundary condition at $x=0$. This shock (at $x=0$) is not physically acceptable and immediately changes to the rarefaction which moves up to the front shock (physically acceptable) moving with the Rankin-Hugoniot speed. The solution endures also after the collision (rarefaction with the front shock). Now we are interested in finding out how many grid points are still sufficient for the numerical approximation. In Fig. 18 we present the solution with the same parameters as in Fig. 17 with only $N=60$ ($M=10$) grid points. The results are nearly the same as the time evolution of concentration with $N=200$ and $M=20$ . This demonstrates the efficiency of our method.

The experiments with higher degeneracy $p<0.85$ are more suitable for our numerical approximations. But when $p\to 1$ ($p=1$ is linear case) troubles arise with our numerical approximations since this is the singular case for our method.

The same mathematical model in 1D was treated numerically in [10] using the operator splitting method, where in the transport part a semi-analytical solution has been used and in the diffusion part the finite volume method has been applied. A very precise numerical solution of this model have been applied in a 2D-problem for the transport of contaminant in a dual-well setting. The jump in the Dirichlet boundary condition is applied there, to create the pulse shape in the injection well and to compute (and also to measure) the response in the extraction well. This gives very important information for the solution of the inverse problems (in the determination of the model parameters). Since it is an ill-posed problem, a very precise numerical solution is required. Therefore our method can be successfully applied there and we will implement it in our forthcoming paper. The method used in [10] is limited to Freundlich and Langmuir isotherms. Our method also works in a much more general case of sorption isotherms, since it does not depend on their geometrical properties. For example, we consider $\Psi(s)=\frac{as^{p}}{1+bs^{p}}$ with $p\in(0,1)$. This isotherm equals to the Freundlich one for $b=0$ and to the Langmuir one for $p=1$ . In this setting our model (37) admits the solution with interface ($p\in(0,1)$). The case $p=1$ is reduced to the Langmuir sorption isotherm which is a nondegenerated problem. This case is a singular one for our method (the speed of interface is $\infty$). However, we take $p=0.95$ and expect that our solution will be close to the one with Langmuir sorption isotherm solved in [10] with $D=0$ (only transport). To compare these results, we take $D=0.005$, $a=1.5$ and $b=1$. We use $N=100$ and $M=20$ in $D_{4}$ discretization. The corresponding results are in Fig. 19, which we can compare (graphically) with the corresponding results in [10] (Fig. 4). Our concentration profiles at time sections 5 and 11 then correspond to the ones for $t=16$ and $t=40$ in [10]. There is a good agreement. This gives us the possibility to obtain a good approximation also in the case of sorption isotherms $\Psi(s),\ |\dot{\Psi)}^{\prime}(0)|<\infty$ representing a nondegenerated problem (without interface). In that case, in the place of $\Psi(s)$ we consider an approximation $\Psi(s^{p})$ with $p\nearrow 1$, which is of course limited by the numerical stability of the ODE solver in our setting. We observe that the $p=0.95$ leads to a good approximation.

In the next experiment we take $b=0.01$ and the remaining parameters are the same as in the previous experiment (with $p=0.95$). The results are shown in Fig. 20 and they represent a good approximation of linear sorption (case $p=1$). Consequently, when a sorption isotherms $\Psi(s)$ doesn’t leads to the appearance of the interface (nondegenerated case), we can use the approximation based on the Langmuir sorption type $\Psi(\frac{s^{p}}{1+bs^{p}})$ with $p\nearrow 1,\ b\searrow 0$ close to the instability region.

Now we compare our results with the corresponding ones in [24] (based on the operator splitting method) with the following parameters: $p=0.75,D=0.1,V=1$ and $a=2$ . We will use the discretization parameters $N=100$ and $M=20$. Time evolutions of the concentrations are presented in Fig 21. The graph corresponding to the $5$-th time section can compared with the one in [24] (Fig. 1 for $t=15$). Then, we use $D=0.001$ and other parameters are the same as in the previous experiment. Then, we can compare the graph of the solution shown in Fig. 22 ($5$-th time section) with the corresponding results in [24] ( Fig. 2 for $t=15$ ). In both cases we obtain a very good agreement.

7.4 Simplified model for turbulent flow in porous media

Consider the following mathematical model

[TABLE]

which is a simplified 1D mathematical model for a turbulent flow in porous media -see [16]. Also here for a special initial profile $u_{0}$, there is an analytical solution given by [16]

[TABLE]

where $(A)_{+}=A$ for $A>0$ otherwise $(A)_{+}=0$. In this case the interface is given by the formula

[TABLE]

The analytical solution satisfies our interface model given by Theorem 1 (similarly as in Section 7.2).

Following our approximation method in Sections 2-4, we successively obtain $\beta=2$ and the ODE-model for the interface evolution ($x\nearrow s(t)$)

[TABLE]

We present our results with discretization parameters $N=60$ and $M=20$ in the discretization strategy $D_{4}$. The time evolution of the numerical solution is in Fig. 23. The analytical solution graphically coincides with the numerical one in the used scaling. The interfaces are presented in Fig. 24 and the errors in Fig. 25.

As we can see the numerical and analytical interfaces cannot be distinguished in the used scaling. The solution is growing since the source term is positive (the term $u^{1/2}$ dominates $u^{3/2}$).

We also present the solution of a foam drainage model (see [15, 26, 27])

[TABLE]

We consider the boundary condition $\partial_{x}u|_{x=0}=0$ (symmetry at $x=0$) and the initial condition $u(x,0)=u_{0}(x)$.

The numerical solution is presented in Fig. 26.

Finally, we present the numerical solution for the mathematical model of viscous liquid advancing over a dry bed. The motion of a thin sheet of viscous liquid over an inclined plate (see [5]) is modelled by

[TABLE]

We consider the boundary condition $\partial_{x}u|_{x=0}=0$ (symmetry at $x=0$) and the initial condition $u(x,0)=u_{0}(x)$.

The numerical solution is presented in Fig. 27.