Numerical approximation of the insitu combustion model using the nonlinear mixed complementarity method

several mathematical models in different disciplines such as engineering, physics, economics and other sciences study partial differential equations of the parabolic type. These models can lead to the problem of mixed complementarity, that is, the case several mathematical models in different disciplines, such as engineering, physics, economics and other sciences study partial differential equations of the parabolic type. These models can lead to the problem of mixed complementarity, that is, the case of the in-situ combustion model, which will be our model in this work. Other applications of complementarity problems are described in [9].

Since the goal is to find an approximation of the analytical solution, we will develop a numerical method that will allow us to achieve our goal. This technique will be applied to the simple in-situ combustion model which will be written as a mixed complementarity problem.

The in-situ combustion model is a particular case of the model treated in [4]. In this case, the model considers the injection of air into a porous medium containing solid fuel and consists of a system of two nonlinear parabolic differential equations.

The contribution of the work is the study of the simple in-situ combustion model and simulations for the proposed model applying the Crank-Nicolson method and the FDA-MNCP algorithm [6].

We will present the results that indicate that the sequence of feasible points generated is contained in a feasible region and we will verify that the directions obtained are feasible and descending for a function associated with the complementarity problem and we will also see the proof of global convergence for the FDA-MNCP following the feat for FDA-NCP[11]. We apply the FDA-MNCP method to the in-situ combustion problem, describe the discretization procedure using the finite difference technique for the mixed complementarity problem associated with the problem, and also present the numerical results and the corresponding error analysis with the comparison with the FDA-NCP method. Finally, we will present some conclusions.

The model studies one-dimensional flows, with a combustion wave in the case when the oxidant (air with oxygen) is injected into a porous medium. Initially the medium contains a fuel that is essentially immobile, does not vaporize and the amount of oxygen is unlimited. This is the case for solid or liquid fuels with low saturations. As in [1], we study the simplified model where:

• $\bullet$

  A small part of the available space is occupied by fuel.

• $\bullet$

  Porosity changes in the reaction are negligible.

• $\bullet$

  Temperature of solid and gas are the same (local thermal equilibrium).

• $\bullet$

  Gas velocity is constant.

• $\bullet$

  Heat loss is negligible.

• $\bullet$

  Pressure variations are small compared to the prevailing pressure.

The model has temporal $t$ and spatial $x$ coordinates that include the heat balance equation, the molar balance equation for immobile fuels, and the ideal gas law.

where $T[K]$ is the temperature, $\rho[\frac{mol}{m^{3}}]$ is the molar density of the gas, and $\rho_{f}[\frac{mol}{m^{3}}]$ is the molar concentration of the immobile fuel. The set of parameters along with their typical values are given in the Table1.

As in [1], if we consider for simplicity $\mu_{f}=\mu_{g}=\mu_{0}=1$ and the amount of oxygen is unlimited, the reaction ratio $W_{r}$ is taken as:

The variables to be found are the temperature $(T)$ and the molar concentration of the immobile fuel $(\rho_{f})$. Since the equations are not dimensioned, we do as [1], to obtain the dimensionless form:

with the dimensionless constants:

Here $P_{e_{T}}$ is the Peclet number for thermal diffusion, $u$ becomes the dimensionless thermal wave velocity, $E$ is a escalated activation energy and $\theta_{0}$ is a scaled reservoir temperature. With a reservoir initial conditions:

and injections condition:

We now describe in detail the finite difference scheme for the in-situ combustion model for the FDA-MNCP method. For this we use the mesh applying the Crank-Nicolson method to approximate the spatial derivatives in each case, i.e.:

Considering the Dirichlet conditions at the point $x_{0}$:

and Neumann conditions at the point $x_{M}$

given in[20], therefore, the value at $x_{0}$ is known at all times but not at $x_{M}$, and

The boundary condition on $x_{M}$ gives:

therefore,

given shaping the FDA-MNCP method:

To obtain the discrete relations of(14), we replace (8-11) in(14) to obtain:

where $\lambda=\frac{h}{h}$ and $\mu=\frac{k}{hr}$. The scheme is valid for $m=1,2,\dots,M$ of the points whose values are not known. At the boundary points we have that for $m=1$, substituting (12) in (16) we obtain:

for $n=M$, we replace (13) in (16) for obtain:

Then (16) is valid for all $m=2,\dots,M-1$ and joining the expressions (17) and (18)we obtain the following inequality in the variable $\theta^{n+1}$

where $LD=B\theta^{n}-\lambda P(\theta^{n},\eta^{n})+2k\Phi(\theta^{n},\eta^{n})+UR$ is known at every instant of time.

Furthermore:

where $A,B\in\mathbb{R}^{M\times M};\ \theta^{n},P^{n},\Phi^{n}\in\mathbb{R}^{M}.$ Similarly, to obtain the discrete form of (15) we replaced (8-11) in (15) to obtain:

The difference scheme is valid for all $m=1,2,\dots,M$ of the points whose values are not known.

At the boundary points we have that for $m=1$, we replace (12) in (24) we obtain:

for $m=M$, we replaced (13) in (24) to obtain:

Therefore, (24) is valid for all $m=2,\dots,M-1$ and joining the expressions (25) and (26) we obtain the following inequality in the variable $\eta^{n+1}$;

where $LDQ=Diag(2)\eta^{n}+k\varphi(\theta^{n},\eta^{n})$ is known at every instant of time.

Furthermore:

Therefore, the discrete form of (14) and (15) is given by (19) and (27).

and must comply with:

Thus, joining (19), (27) (30) and (32) form a Mixed Complementarity Problem, which can be solved by the FDA-MNCP algorithm

The numerical results obtained from the implementation of the Algorithm in Matlab are shown in the following Section.

We will now present the numerical results of the simulations performed in Matlab for the FDA-MNCP method. For this simulation, we considered $[x_{0},x_{M}]=[0,\ 0.05]$ and $[t_{0},t_{N}]=[0,1]$ as the space and time intervals respectively. We keep the numbers of subintervals in time constant and equal to: $N=10^{5}$, that is, $k=\Delta t=10^{-5}$ while the number of subintervals in space will be equal to $M=50,100,200,400$. For the FDA-MNCP method we consider an error tolerance of $10^{-8}$.

The values of the dimensionless parameters in (7) are:

with the previous input data we obtain Figures (1), (2), (3), (4), which show the results obtained by algorithm 3.1 and algorithm 5 of [1] for the FDA-MNCP and FDA-NCP [11] methods, respectively.

As we can see in Figures 1, 2, 3 and 4 obtained previously, they show us that the results obtained by the FDA-MNCP methods and FDA-NCP, coincide very well, as can be seen in the times indicated in the figures presented previously. In Figure 5, the differences between the methods can be observed.

we can see that the difference between the solutions of $\theta$ and $\eta$ are very small as we increase the number of points, between the FDA-MNCP and FDA-NCP methods[11].

Below we show four tables where we compare the computational process time for $M=50,\ 100,\ 20,\ 400$ of the FDA-MNCP method and the FDA-NCP method studied in [1].

As we can see in Table 2, the computational process time used by the FDA-MNCP method doubles the time of the FDA-NCP method for a partition of the $x$-axis into fifty points; the number of iterations of the FDA-MNCP method for the indicated times is also greater than that of the FDA-NCP method. In addition, the third column of each table shows the value of $t^{k}$ for each of the methods.

In Table 3, we can also observe that the computational process time used by the FDA-MNCP method is slightly greater than the time of the FDA-NCP method as the number of iterations.

In Table 4, we see that when we increase the number of partitions to $200$, the computational process time of the FDA-MNCP method starts to be shorter than the time of the 1FDA-NCP method, although the number of iterations is still greater.

With the last Table 5, we can observe that by considerably increasing the number of partitions, the computational process time of the FDA-MNCP method is much lower than that of the FDA-NCP[11] method, because for the calculation of $\nabla S$ in the FDA-MNCP method, it makes half the calculations than for $\nabla F$ in the FDA-NCP[11] method, although the number of iterations is greater. We can also observe that the number of iterations varies little in each of the tables presented previously.

In this section we will perform a numerical study of the relative error for the FDA-MNCP method, for each instant of time and $E_{\Delta x}$, where $\Delta x\ =\ h,\ \ h/2,\ \ h/4$. The length of the subinterval in time will be constant and equal to $\Delta t\ =\ k\ =\ 10^{-5}$  e  $h\ =\ \frac{1}{50}$, and compared with the FDA-NCP method, [1]. Tables 6 and 8 show the results for the Relative Errors of the FDA-MNCP method, while Tables 7 and 9 show the Relative Error for the FDA-NCP method and are represented in Figures 6 and 7.

As can be seen in Figures 6 and 7, the relative errors corresponding to the FDA-MNCP and FDA-NCP methods, respectively, are very similar, as can be seen in the following tables, which show the relative errors for $\theta$ and $\eta$ usando $h=1/50$.

In Tables 6 and 7 we see that the relative errors $E_{h},\ E_{\frac{h}{2}},\ E_{\frac{h}{4}}\$ for $\ \theta\$ of both the FDA-MNCP method and the FDA-NCP method are very similar since they are the same in the first three decimal places.

Similarly, in Tables 8 and 9 we $E_{h},\ E_{\frac{h}{2}},\ E_{\frac{h}{4}}\$ for $\eta$ of both the FDA-MNCP method and the FDA-NCP method are very similar, since the first three decimal places are the same

• •

  We propose the solution of the simple in-situ combustion model using a numerical method based on an implicit finite difference scheme and a nonlinear mixed complementarity algorithm, which can be applied to parabolic problems and can be rewritten in the form of a mixed complementarity problem. This method is applied to the system(1).

• •

  In this work, we have seen that the feasible interior point algorithm FDA-MNCP, as a good technique to numerically solve mixed complementarity problems. Theoretical results of the algorithm ensure the global convergence.

• •

  We solve the in-situ combustion model(1) using the nonlinear mixed complementarity algorithm and comparing it to the FDA-NCP method, with the two solutions being very close, as can be seen in Figures 1, 2, 3 and 4. This suggests that it is possible to apply this method to parabolic and hyperbolic problems, which can be written as a mixed complementarity problem. The FDA-MNCP method shows the advantage of being able to be used with more points in the discretization of the space, as can be seen in Tables 4 and 5, in addition to being faster in computational time than the FDA-NCP method when we increase the discretization in the space.

• •

  Regarding the relative errors, Tables 6 and 8 show good evidence of convergence of the mixed complementarity algorithm, which is observed in Figure 5, which indicates a decrease in linear growth when we refine the mesh.