Numerical Simulation of Production Localization Problem with Risks of Technological Disaster

This article develops numerical algorithms for solving inverse problems of recovering diffusion coefficient from additional information about the solution of reaction-diffusion equation. The results of computational experiments are presented and discussed. Their contribution to the methods of mathematical modeling of technological disasters and the solution of problems of peoples safety ensuring in case of living near potentially dangerous objects is estimated.


Introduction
This article continues a small cycle of papers [1][2][3][4] devoted to the theoretical and numerical study of boundary and extremun problems for the reaction-diffusion equation and their applications. The reaction-diffusion model describes the process of the pollutant propagation, in particular, the release and spread of substances which are hazardous to human health after a techogenic accident. The articles [1,2] consider the application of these models while solving problems of the location of enterprises, where technological disasters could be possible. It was assumed that the accident was inevitable and the volumes of pollutant emissions were known. Using mass transfer models, we establish the possible concentration of pollutants in the domain of interest. The latter is determinative when choosing both the location of the enterprise itself and its infrastructure. Let us note that the mass transfer model contains a number of physical parameters that may not be known. In this case, inverse problems arise, which consist in finding the indicated parameters from additional information about the solution of the boundary value problem [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].
Unlike [1,2], this article focuses on numerical algorithms for solving inverse problem of diffusion coefficient's recovering from the pollutant measured in a certain subdomain. This problem is very important for applications. Assuming that the diffusion coefficient depends only on a spatial variable or is constant, we can find it by measuring the concentration of a (harmless)substance in a measurable domain. Obviously, knowing the parameters of the model, we can easily solve the direct problem. Finally, we note that in this paper we develop numerical algorithms based on optimality systems for inverse coefficient problems obtained in [1,2].

Boundary value and extremum problems
We study the following mass (pollutant) transfer model: (1) considering in a bounded domain Ω of the space d R , d = 2,3 with Lipschitz boundary Γ. Here λ≡λ(x)>0 is the variable diffusion coefficient, k≡k(x)≥0 is quantity characterizing disintegration of pollutant by chemical reactions, f≡f(x) is the volume source density, ψ=ψ(x) is a function given on boundary Γ.
As we can use this information, for example, concentration d(x) is measured in some subdomain QΩ. Let coefficients λ is unknown function and we must determine this function together with the solution  of problem (1).
Our inverse extremum problem consists of finding two functions (φ, λ) such that here μ0 ≥ 0, μ1 ≥ 0 are positive constants. Description of the other symbols used can be found, for example, in [17,Chapter 3].
Using the mathematical apparatus of the book [17] we obtained the optimality system. This system has the meaning of the necessary extremum condition of the first-order. The optimality system is essentially used to prove the uniqueness and stability of the solution of the extremum problem and for creation of numerical algorithms. This system consists of the direct problem, adjoint problem and variational inequality which has the meaning of the minimum principle. The optimality system has the following form: . 0

Results of numerical experiments
Similar to papers [17][18], algorithm of numerical solution of inverse extreme problem (2) based on system of optimality (3), (4) and (5) was built. Discretization of boundary problems was carried out using the finite differences method, and the solution to the optimality system was based on using the Newton method [20].The algorithm was tested using the Scilab software package [21].
The initial data was set as follows. The domain Ω was selected as a unit square (see Fig. 1).The source of pollution was described by the formula At the first stage of the study, the effect of the regularization parameter μ1 on solving the inverse extreme problem was studied. Figure 2 and figure 3 show the behavior of the relative errorsE0 and E1 for the two initial approximation choices. Figure 2    The second stage of research was to determine the number of iterations in the Newton method that will be required to obtain a solution with a certain accuracy. Figure 4 and figure 6 show graphs of the dependencies of errors E0 and E1 on the number of iterations for each of the initial approximations λ1 and λ2.The graphs show that 5-6 iterations are sufficient to restore the values of the concentration  function with an accuracy of 10 -6 , while the coefficient λ is restored with an accuracy of 10 -4 . Figure 5 and Figure 7 show the process of restoring the coefficient λ in the form of slices of the surfaces λd ("Toch") and λ ("Iter = [1:5]") at y =0.5. Graphs are shown for 5 iterations, since for the following iterations the differences in the graphs are invisible.    It is well known that Newton's method requires a good initial approximation. And many researchers wonder what initial approximation to choose when solving a particular problem. As mentioned above, as an initial approximation, the functions λ1 and λ2 were chosen. The function λ1 is the simplest choice of the initial approximation and at the same time repeats the form of the sought function λ, while the form of the function λ2 is significantly different from the sought function, but at the same time makes it possible to find the sought function λ with an accuracy of µ1=10 -8 . To obtain a group of functions that are initial approximations for task (2), the parameter b was used, which is part of the selected functions λ1 and λ2. Figure 8 (for λ1) and Figure 9 (for λ2) show the value ranges of parameter b, using the relative error graphs E0 and E1.  Figure 8.  The function responsible for the concentration of the substance is restored with accuracy 10 -6 , and the diffusion coefficient with accuracy 10 -3 (see Figure 11 and Figure 13).The maximum error value at a point close to (0,0) is associated with an inequality (5) in which there is a term    . At the specified point, it quickly tends to 0 and, unfortunately, does not allow you to qualitatively restore the values of the functions you are looking for at this point.     Based on the results of the studies, it can be concluded that the developed algorithm for solving the inverse extreme problem for the diffusion-reaction model (2) allows you to solve quickly these problems and restore unknown parameters with sufficiently high accuracy. By changing the parameters of the algorithm, such as the regularization parameter or the number of iterations, it is possible to obtain the required accuracy in solving inverse extreme problems. Initial approximations can be functions that both repeat the geometry of the desired parameters and have different geometry.

Conclusion
In this work, the main attention is paid to the problem of restoring the diffusion coefficient depending on the spatial variable. For this inverse problem the efficient numerical algorithm was developed and implemented, based on the optimality system. The algorithm was implemented in a package Scilab. The study of various parameters of the developed algorithm showed its effectiveness in finding the diffusion coefficient and concentration of matter. The obtained algorithm for finding unknown parameters for the diffusion-reaction model will allow modeling the process of distribution of pollutants resulting from a man-made disaster. Analysis of simulation results will allow to identify areas exposed to high level of contamination and make them non-residential at the stage of production design.