The work is dedicated to mathematical modelling of random diffusion flows of admixture particles in a two-phase stratified strip with stochastic disposition of phases and random thickness of inclusion-layers. The study of such models are especially important during the creation of composite layered materials, in the research of the transmission properties of filters, and in the prediction of the spread of pollutants in the environment. Within the model we consider one case of uniform distribution of coordinates of upper boundaries of the layers of which the body is made up and two more cases, i.e., of uniform and triangular distributions of the inclusion thickness. The initial-boundary value problems of diffusion are formulated for flux functions; the boundary conditions at one of the body’s surfaces are set for flux and, at the other boundary, the conditions are given for admixture concentration; the initial condition being concerned with zero and non-zero constant initial concentrations. An equivalent integro-differential equation is constructed. Its solution is found in terms of Neumann series. For the first time it was obtained calculation formulae for diffusion flux averaged over the ensemble of phase configurations and over the inclusion thickness. It allowed to investigate the dependence of averaged diffusion fluxes on the medium’s characteristics on the basis of the developed software. The simulation of averaged fluxes of admixture in multilayered $Fe-Cu$ and $\alpha Fe-Ni$ materials is made. Comparative analysis of solutions, depending on the stage of averaging procedure over thickness, is carried out. It is shown that for some values of parameters the stage of averaging procedure over thickness has almost no effect on the diffusion flow value.

We propose and justify a numerical method for computing the double integral with variable upper limits that leads to the variableness of the region of integration. Imposition of simple variables as functions for upper limits provides the form of triangles of integration region and variable in the external limit of integral leads to a continuous set of similar triangles. A variable grid is overlaid on the integration region. We consider three cases of changes of the grid for the division of the integration region into elementary volumes. The first is only the size of the imposed grid changes with the change of variable of the external upper limit. The second case is the number of division elements changes with the change of the external upper limit variable. In the third case, the grid size and the number of division elements change after fixing their multiplication. In these cases, the formulas for computing double integrals are obtained based on the application of cubatures in the internal region of integration and performing triangulation division along the variable boundary. The error of the method is determined by expanding the double integral into the Taylor series using Barrow’s theorem. Test of efficiency and reliability of the obtained formulas of the numerical method for three cases of ways of the division of integration region is carried out on examples of the double integration of sufficiently simple functions. Analysis of the obtained results shows that the smallest absolute and relative errors are obtained in the case of an increase of the number of division elements changes when the increase of variable of the external upper limit and the grid size is fixed.