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

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

$$\partial_{t}u=\partial_{x}^{2}u^{n}+b_{0}\partial_{x}u^{\gamma}+c_{0}u^{m}$$

(1)

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

$$a)\quad u=\phi(t),\ \mbox{or}\quad b)\quad-\partial_{x}u^{n}+b_{0}u^{\gamma}=q(t).$$

(2)

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

$$u(x,0)=u_{0}(x)\geq 0,$$

(3)

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

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

(4)

is of the form

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

(5)

where

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

and the corresponding initial condition is

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

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

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

(6)

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

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

(7)

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.

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

$$\bar{u}_{t}=\bar{{\cal L}}(\bar{u})-\frac{\dot{s}}{s}y\partial_{y}\bar{u},\ y\in(0,1)$$

(8)

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))

$$\dot{s}={\cal B}(\partial_{y}\bar{u}(1,t)),\ \mbox{resp.}\ \dot{s}={\cal B}(\partial_{y}\bar{u}(1,t),\partial_{y}^{2}\bar{u}(1,t)),$$

(9)

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.

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

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

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

$$-\sigma f^{\prime}=nf^{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,\quad\xi\in(-\infty,\infty).$$

The function $\theta$ defined as

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

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

$$\theta(s)=\sigma s+b_{0}s^{\gamma}-nc_{0}\int_{0}^{s}\frac{r^{m+n-1}}{\theta(r)}dr.$$

(10)

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

$$\theta(s)\approx\theta_{0}s^{q}\quad\mbox{as}\ s\searrow 0$$

(11)

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

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

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

$$\theta(s)=\theta_{0}s^{q}+o(s^{q}),\quad\mbox{with}\ s^{-q}o(s^{q})\to 0,\ \mbox{for}\ s\searrow 0.$$

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)

$$f(\xi)\approx d^{\beta}(\xi_{0}-\xi)^{\beta}\quad\mbox{for}\ \xi\nearrow\xi_{0}.$$

(12)

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

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

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

$$f(\xi)\approx[\frac{\theta_{0}(n-q)}{n}]^{\beta}{(\xi}_{0}-\xi)^{\beta},\quad\mbox{where}\ \beta=\frac{1}{n-q},d=\frac{\theta_{0}(n-q)}{n}$$

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

$$f(\xi)\approx d^{\beta}(\xi_{0}-\xi)^{\beta},\ f^{\prime}(\xi)\approx\beta d^{\beta}(\xi_{0}-\xi)^{\beta-1}\quad\mbox{for}\ \xi\nearrow\xi_{0},$$

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

(13)

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

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

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

(14)

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

$$\dot{s}(t)=-\frac{\partial_{t}u}{\partial_{x}u}\ \mbox{in the sense}\ x\nearrow s(t),$$

(15)

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

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

(16)

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

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

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

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

$$|\frac{w(x,t)-w(s(t),t)}{s(t)-x}-d|\to 0,\ \mbox{as}\ x\nearrow s(t).$$

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

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

$$\partial_{t}w=n(n\beta-1)w^{(n-1)\beta}\partial_{x}^{2}w+nw^{(n-1)\beta-1}(\partial_{x}w)^{2}+b_{0}\gamma w^{(\gamma-1)\beta}\partial_{x}w+\frac{c_{0}}{\beta}w^{(m-1)\beta+1}.$$

(17)

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

$$\partial_{t}w=n(n\beta-1)\frac{w^{(n-1)\beta}}{s^{2}}\partial_{y}^{2}w+n\frac{w^{(n-1)\beta-1}}{s^{2}}(\partial_{y}w)^{2}+b_{0}\gamma\frac{w^{(\gamma-1)\beta}}{s}\partial_{y}w+$$

$$+\frac{c_{0}}{\beta}w^{(m-1)\beta+1}+y\frac{\dot{s}}{s}\partial_{y}w.$$

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

$$\partial_{t}w=n(n\beta-1)\frac{w^{(n-1)\beta}}{s^{2}}\partial_{y}^{2}w+n\frac{w^{(n-1)\beta-1}}{s^{2}}(\partial_{y}w)^{2}+b_{0}\gamma\frac{w^{(\gamma-1)\beta}}{s}\partial_{y}w+$$

$$+\frac{c_{0}}{\beta}w^{(m-1)\beta+1}+y\frac{Q(w)}{s}\partial_{y}w,$$

(18)

$$\dot{s}=Q(w)$$

(19)

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

$$\partial_{y}w(1,t)=\frac{d\ell_{B}}{dy}|_{y=1},$$

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

$$\dot{C}_{i}=n(n\beta-1)\frac{C_{i}^{(n-1)\beta}}{s^{2}}\frac{d^{2}\ell_{i}}{dy^{2}}|_{y=y_{i}}+n\frac{C_{i}^{(n-1)\beta-1}}{s^{2}}(\frac{d\ell_{i}}{dy}|_{y=y_{i}})^{2}+b_{0}\gamma\frac{C_{i}^{(\gamma-1)\beta}}{s}\frac{d\ell_{i}}{dy}|_{y=y_{i}}+$$

(20)

$$+\frac{c_{0}}{\beta}C_{i}^{(m-1)\beta+1}+y_{i}\frac{Q(C)}{s}\frac{d\ell_{i}}{dy}|_{y=y_{i}},$$

$$\dot{s}=Q(C),\quad\mbox{where}\ \partial_{y}w(1,t)\leftrightarrow\frac{d\ell_{B}}{dy}|_{y=1}$$

(21)

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

$$C_{-1}=C_{1}+s(-b_{0}C_{0}^{\gamma\beta}+q)\frac{2\alpha_{1}}{n\beta C_{0}^{n\beta-1}}.$$

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.

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

$$a(t,u)=g(t,u)u^{n},\ b(x,t,u)=b_{0}(x,t)u^{\gamma},\ c(x,t,u)=c_{0}(x,t)p(u)u^{m}$$

(22)

$$\mbox{with}\quad g(t,0)\neq 0,\ p(0)\neq 0\quad\mbox{and}\ g\in C^{2}.$$

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

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

(23)

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

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

$$\partial_{t}u=D\partial_{x}^{2}u-H(u),\ x\in(0,\infty),t>0$$

(24)

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

$$H(z)=\left\{\begin{array}[a]{c c}1&\ \mbox{for}\ z>0\\ \\ 0&\ \mbox{for}\ z=0\end{array}\right.$$

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

$$a=a_{0}(x,t)u,\ b=b_{0}(x,t)u,\ c=c_{0}(x,t)H(u)u^{m}$$

(25)

$$\mbox{with}\quad a_{0}(x,t)\geq\delta>0,\ c_{0}(x,t)\leq-\delta,\ 0\leq m<1.$$

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

$$\dot{s}(t)=-\frac{\partial_{x}^{2}a(x,t,u)+\partial_{x}b(x,t,u)+c(x,t,u)}{\partial_{x}u},$$

(26)

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

$$\partial_{t}w=a_{0}[\partial_{x}^{2}w+(\alpha-1)\frac{\partial_{x}w)^{2}}{w}]+[b_{0}+2\partial_{x}a_{0}]\partial_{x}w+\frac{1}{\alpha}\frac{w^{(m-1)\alpha+2}c_{0}}{w}+\frac{1}{\alpha}(\partial_{x}^{2}a_{0}+\partial_{x}b_{0})w$$

(27)

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

$$(\alpha-1)a_{0}(\partial_{x}w)^{2}+\frac{c_{0}}{\alpha}w^{(m-1)\alpha+2}\to 0\quad\mbox{for}\ x\to s(t).$$

(28)

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

$$\partial_{x}w=-\sqrt{-\frac{c_{0}}{\alpha(\alpha-1)a_{0}}w^{(m-1)\alpha+2}}\quad\mbox{for}\ x\nearrow s(t).$$

(29)

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

$$\dot{s}(t)=-\frac{\partial_{t}u}{\partial_{x}u}=-\frac{\partial_{t}w}{\partial_{x}w}$$

(30)

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

$$(m-1)\alpha+2=0\quad\Rightarrow\ \alpha=\frac{2}{1-m}>1.$$

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

$$Z_{w}=\lim_{x\nearrow s(t)}\left(\frac{(\alpha-1)a_{0}(\partial_{x}w)^{2}+\frac{1}{\alpha}c_{0}}{w}\right).$$

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

$$Z_{w}=(\alpha-1)\partial_{x}a_{0}(s(t),t)\partial_{x}w+2(\alpha-1)a_{0}(s(t),t)\partial_{x}^{2}w+\frac{1}{\alpha}\frac{\partial_{x}c_{0}(s(t),t)}{\partial_{x}w}=$$

(31)

$$2(\alpha-1)a_{0}(s(t),t)\partial_{x}^{2}w-$$

$$(\alpha-1)\partial_{x}a_{0}(s(t),t)\sqrt{-\frac{c_{0}(s(t),t)}{\alpha(\alpha-1)a_{0}(s(t),t)}}-\frac{1}{\alpha}\partial_{x}c_{0}(s(t),t)\sqrt{-\frac{\alpha(\alpha-1)a_{0}(s(t),t)}{c_{0}(s(t),t)}}.$$

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

$$\dot{s}(t)=(2\alpha-1)a_{0}(s(t),t)\sqrt{-\frac{\alpha(\alpha-1)a_{0}(s(t),t)}{c_{0}(s(t),t)}}\partial_{x}^{2}w(s(t),t)-b_{0}(s(t),t)+$$

(32)

$$2\partial_{x}a_{0}(s(t),t)+(\alpha-1)\partial_{x}c_{0}(s(t),t)\frac{a_{0}(s(t),t)}{c_{0}(s(t),t)},\ \mbox{where}\ \alpha=\frac{2}{1-m}\ \mbox{and}\ u=w^{\alpha}.$$

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

$$(y_{N-3},C_{N-3}),(y_{N-2},C_{N-2}),(y_{N-1},C_{N-1}),(1,0).$$

has to be used.