*Article* **Numerical Study of Soil-Thawing Effect of Composite Piles Using GMsFEM**

**Petr V. Sivtsev 1,\*, Piotr Smarzewski <sup>2</sup> and Sergey P. Stepanov <sup>1</sup>**


**Abstract:** During construction works, it is advisable to prevent strong thawing and an increase in the moisture content of the foundations of engineering structures in the summer. Since the density of water and ice differ, due to the difference bulging of the foundation sections can occur when it freezes back in winter. In this work, the effect of fiber-reinforced piles on the thermal field of the surrounding soil is investigated numerically; that is, the study of the influence of aggregates with high and low thermal-physical properties on the temperature of frozen soils is conducted. Basalt and steel fiber reinforcement are compared. The difficulty of this work is that the inclusions inside piles are too small compared to the pile itself. Therefore, to solve the Stefan problem, a generalized multiscale finite element method (GMsFEM) was used. In the GMsFEM, the usual conforming partition of the domain into a coarse grid was used. It allowed reducing problem size and, consequently, accelerating the calculations. Results of the multiscale solution were compared with fine-scale solution, the accuracy of GMsFEM was investigated, and the optimal solution parameters were defined. Therefore, GMsFEM was shown to be well suited for the designated task. Collation of basalt and steel fiber reinforcement showed a beneficial effect of high thermal conductive material inclusion on freezing of piles in winter.

**Keywords:** Stefan problem; multiscale; generalized multiscale finite element method; composite pile; thermal conduction

### **1. Introduction**

Thermal calculations are important in the construction of engineering geotechnical structures and buildings in the permafrost zone. The temperature regime (a set of sequential temperature fields in the soil mass corresponding to any given points in time from the beginning of the calculation) is calculated as the forecast of the thermal effects on the upper and lower boundaries of the structure foundation set for the entire calculation period [1]. The most characteristic feature of these processes is the previously unknown ("free") boundaries between various thawed and frozen states of the soil. Due to this feature, their mathematical models are non-linear and difficult to analyze. The results of calculations of climatic phenomena, which are critical for geotechnical structures on permafrost foundations (in the permafrost zone), show a continuous repetition of the same constant seasonal cycle. This determines periodic harmonic oscillations of all quantities that determine the state of the foundation of a geotechnical structure, including temperature. However, cyclic harmonic changes in the foundation state caused by constantly repeating influences also occur during other processes. The course of temperature change during one period can have a different character. For example, the temperature can change abruptly, continuously increase or decrease, etc.

During construction works, it is advisable to prevent strong thawing and an increase in the moisture content of the foundations of engineering structures in the summer. Since

**Citation:** Sivtsev, P.V.; Smarzewski, P.; Stepanov, S.P. Numerical Study of Soil-Thawing Effect of Composite Piles Using GMsFEM. *J. Compos. Sci.* **2021**, *5*, 167. https://doi.org/ 10.3390/jcs5070167

Academic Editor: Stelios K. Georgantzinos

Received: 22 May 2021 Accepted: 18 June 2021 Published: 28 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the density of water and ice differ, due to the difference bulging of the foundation sections can occur when it freezes back in winter. Repetition of this phenomenon for many years causes the destruction of engineering structures and buildings. In the regions of the Far North, it leads to catastrophic emergencies.

In this work, the effect of fiber-reinforced piles on the thermal field of the surrounding soil is numerically investigated [2], that is, the study of the influence of aggregates with high and low thermal-physical properties on the temperature of frozen soils is conducted [3,4].

Basalt fibers are produced from molten basalt rock. They have very low thermal conductivity, good strength properties and, at the same time, are relatively cheap. Due to these characteristics, they are used in concrete [5–7]. However, in addition to the chemical and mechanical properties of basalt fibers, their cost varies, depending on the type and quality of the raw material and the production process of these fibers [8]. Nevertheless, the characteristics mentioned above and an environmentally friendly manufacturing process [9] might determine their application in high performance concrete structures instead of the most commonly used steel and polypropylene fibers. The length-to-diameter aspect for basalt fibers is about 1000.

As for steel fiber hybridization, industrial or recycled steel fibers can be used [10–13]. According to results available in the scientific literature [14,15], they show similar mechanical responses, both in terms of tensile strength and matrix-to-fiber bond. Recycled steel fibers are derived from waste tires. Their geometrical characterization can be highly variable: they are generally characterized by a nominal diameter ranging between 0.1 and 2 mm with a corresponding average aspect ratio (i.e., length-to-diameter ratio) ranging between 20 and 150. These variations mainly depend on both the original source (i.e., tires typology) and recycling processes.

The volume concentration of fibers is equal to 10%. In our article, we use the lengthto-diameter ratio for both types of fibers equal to 32 for structured and 10 for random distribution of fibers for calculation convenience. Fibers' mechanical properties are presented in Table 1.


**Table 1.** Mechanical properties of fibers.

One of the topical problems of mathematical physics is the Stefan problem, which is an initial–boundary value problem for a parabolic differential equation with discontinuous coefficients, and serves as a mathematical model of the change in the phase state of a substance with unknown interfaces [16,17]. The phase transition boundary is represented as a smearing zone when the phase transition occurs in a given temperature range [16,18].

The difficulty of this work lies in the pile inclusions that are too small compared to the solution area. For such geometrically complex problems, reduction techniques such as multiscale methods are preferred [19,20]. Multiscale methods are becoming very popular at this time [21–23]. Paper [24] develops a multiscale eigenvalue method based on multiscale substructure technology for multiscale analysis of periodic composite structures. This provides a comparative study from a user's viewpoint for multiscale methods, including the mathematical homogenization method, heterogeneous multiscale method, and multiscale finite element method. The article [25] studies the multiscale finite element method for solving elliptic problems arising from composite materials. The influence of production porosity, diameter [26], and density of structures on the thermal conductivity of a composite are considered using multiscale finite element modeling [27]. In [28,29], mathematical modeling is applied for the thermodynamic analysis and design of composite structures.

This article is organized as follows. Section 2 presents the formulation of a mathematical model that describes the dynamics of the temperature distribution with the phase transitions taken into account. In the third, finite element approximation is given. In the fourth, a generalized multiscale finite element method (GMsFEM) is considered. Section 5 presents the numerical results in a two-dimensional setting.

### **2. Mathematical Model**

Let us consider a mathematical model describing the dynamics of the temperature distribution, taking into account the phase transitions of pore moisture into ice and back, at a certain given temperature of the phase transition *<sup>T</sup>*<sup>∗</sup> in area <sup>Ω</sup> <sup>=</sup> <sup>Ω</sup><sup>−</sup> <sup>∪</sup> <sup>Ω</sup>+. Where <sup>Ω</sup>+(*t*) <sup>=</sup> {*x*|*<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>, *<sup>T</sup>*(*x*, *<sup>t</sup>*) <sup>&</sup>gt;*T*} is the area occupied by the thawed soil, where the temperature exceeds the phase transition temperature and Ω−(*t*) = {*x*|*x* ∈ Ω, *T*(*x*, *t*) < *T*} is an area occupied by the frozen ground. A phase transition occurs at the interface between thawed and frozen soils *S* = *S*(*t*) (Figure 1).

**Figure 1.** Phase transition.

To simulate heat transfer processes with phase transitions, the Stefan model is used. It describes thermal processes accompanying phase transformations of the medium with absorption and release of latent heat:

$$\frac{1}{2}\left(c\rho(T) + m\rho\_l L\delta(T - T^\*)\right)\frac{\partial T}{\partial t} - \text{div}(\lambda(T)T) = f,\text{ x} \in \Omega,\ t \in (0, t\_{\text{max}}],\tag{1}$$

where *L* is the specific heat of phase transition, m is a porosity, *δ*(*T* − *T*∗) is Dirac delta function. For the coefficients of the equation, there are following relations:

$$\mathcal{c}\rho(T) = \begin{cases} \mathcal{c}^-\rho^-\prime \; , \; T < T^\*\prime \; , \\\mathcal{c}^+\rho^+\prime \; , \; T \ge T^\*\prime \; , \end{cases} \\ \lambda(T) = \begin{cases} \begin{array}{c} \lambda^-\prime \; , \; T < T^\*\prime \\\ \lambda^+\prime \; , \; T \ge T^\*\prime \end{array}$$

where *c*+, *ρ*+, *λ*<sup>+</sup> and *c*−, *ρ*−, *λ*<sup>−</sup> are the specific heat, density, and thermal conductivity of thawed and frozen soil, respectively.

Since the process of heat propagation is considered in a saturated porous medium, there are following thermal-physical characteristics:

$$c^-\rho^- = (1-m)c\_{\mathfrak{s}\mathfrak{c}}\rho\_{\mathfrak{s}\mathfrak{c}} + mc\_i\rho\_{i\mathfrak{c}}\ c^+\rho^+ = (1-m)c\_{\mathfrak{s}\mathfrak{c}}\rho\_{\mathfrak{s}\mathfrak{c}} + mc\_w\rho\_{w\mathfrak{s}\mathfrak{c}}$$

where *m* is a porosity. The subscripts *sc*, *w*, *i* denote the skeleton of the porous medium, water, and ice, respectively. For the coefficients of the thermal conductivity in the thawed and frozen zones, there are similar relations:

$$
\lambda^- = (1 - m)\lambda\_{\mathfrak{sl}} + m\lambda\_{i\mathfrak{i}} \\
\lambda^+ = (1 - m)\lambda\_{\mathfrak{sl}} + m\lambda\_{\mathfrak{w}\mathfrak{i}}.
$$

In practice, phase transformations do not take place instantly and can proceed in a small temperature range [*T* − Δ, *T* + Δ]. The discontinuous coefficients of Equation (1) are replaced by sufficiently smooth functions of temperature:

$$(c\rho)\_{\Lambda}(T) = c^{-}\rho^{-} + m(c^{+}\rho^{+} - c^{-}\rho^{-})\frac{1}{2}\left(1 + \text{erf}\left(\frac{T - T^{\*}}{\sqrt{2\Delta^{2}}}\right)\right) + \frac{1}{\sqrt{2\pi\Delta^{2}}}e^{-\frac{(T - T^{\*})^{2}}{2\Delta^{2}}},$$

$$\lambda\_{\Lambda}(T) = \lambda^{-} + m(\lambda^{+} - \lambda^{-})\frac{1}{2}\left(1 + \text{erf}\left(\frac{T - T^{\*}}{\sqrt{2\Delta^{2}}}\right)\right).$$

Thus, the following equation for the temperature is obtained:

$$(\varepsilon \rho)\_{\Lambda}(T) \frac{\partial T}{\partial t} - \text{div}(\lambda\_{\Lambda}(T)T) = f, \text{ x} \in \Omega, \ t \in (0, t\_{\text{max}}], \tag{2}$$

Thus, Equation (2) is a multidimensional quasilinear parabolic equation with smooth coefficients.

For piles and fibers the following is true: (*cρ*)<sup>Δ</sup> = (*cρ*)*p*, *<sup>f</sup>* and *λ*<sup>Δ</sup> = *λ*p,f, where *cp*, *<sup>f</sup>* , *ρp*, *<sup>f</sup>* , *λ*p,f are specific heat, density, and thermal conductivity of piles and fibers. That is, there is no phase transition in the pile and fiber areas, respectively.

Equation (2) is supplemented with the initial and boundary conditions

$$T(\mathbf{x},0) = T\_{0\prime}$$

$$-\lambda \frac{\partial T}{\partial \eta} = \mathfrak{a}(T - T\_{air}), \ \mathbf{x} \in \Gamma\_t,$$

$$-\lambda \frac{\partial T}{\partial \eta} = 0, \ \mathbf{x} \in \frac{\Gamma}{\Gamma\_t}.$$

Here *η*—normal vector.

### **3. Fine-Scale Approximation**

The quasilinear parabolic Equation (2) with the corresponding boundary and initial conditions is approximated using the finite element method in combination with a purely implicit linearized finite difference approximation in time. This means that during discretization, the coefficients depending on the desired function are taken from the previous time layer. Let us write down the variational problem statement for each time layer: find *<sup>T</sup>* ∈ *<sup>H</sup>*<sup>1</sup> such that:

$$a\left(T^{n+1}, v\right) = L(v), \forall \in H\_0^1.$$

where

$$\begin{array}{rcl} a\left(T^{n+1}, \upsilon\right) &=& \frac{1}{\tau} \int\_{\Omega} (c\rho)\_{\sigma}(T^{n})T^{n+1}w dx + \int\_{\Omega} \lambda\_{\sigma}(T^{n})T^{n+1} \cdot w dx + \int\_{\Gamma t} aT^{n+1} ds, \\ L(\upsilon) &=& \frac{1}{\tau} \int\_{\Omega} (c\rho)\_{\sigma}T^{n}w dx + \int\_{\Gamma t} aT\_{air} ds. \end{array}$$

To solve the problem numerically, it is necessary to pass from a continuous variational problem to a discrete one. For this, finite-dimensional spaces *<sup>V</sup>*<sup>h</sup> <sup>∈</sup> *<sup>H</sup>*1, *<sup>V</sup>*<sup>ˆ</sup> *<sup>h</sup>* ∈ *<sup>H</sup>*<sup>1</sup> <sup>0</sup> are introduced and the following problem is defined on them: find *Th* ∈ *Vh* such that:

$$a\left(T\_h^{\mathfrak{u}+1}, \upsilon\right) = L(\upsilon), \,\forall \,\in \,\hat{\mathbb{V}}\_{h\nu}$$

where

$$\begin{array}{rcl} a\left(T^{n+1}, \upsilon\right) &=& \frac{1}{\tau} \int\_{\Omega} (c\rho)\_{\sigma} (T^{n}) T^{n+1} w\_{\hbar} d\mathbf{x} + \int\_{\Omega} \lambda\_{\sigma} (T^{n}) T^{n+1}\_{\hbar} \cdot w\_{\hbar} d\mathbf{x} + \int\_{\Gamma t} a T^{n+1} d\mathbf{s}, \\\\ L(\upsilon) &=& \frac{1}{\tau} \int\_{\Omega} (c\rho)\_{\sigma} T^{n}\_{\hbar} w\_{\hbar} d\mathbf{x} + \int\_{\Gamma t} a T\_{\mathrm{air}} d\mathbf{s}. \end{array}$$

### **4. Generalized Multiscale Finite Element Method (GMsFEM)**

Construction local reduction of a model on the snapshot space is described by solving several local spectral problems using GMsFEM. In the GMsFEM the usual conforming partition of domain into finite elements is used. This partition is called the coarse grid T*H*. The coarse grid is partitioned into coarse-grid blocks. That is, T*<sup>H</sup>* is a coarse grid in domain <sup>Ω</sup>, such that <sup>T</sup>*<sup>H</sup>* <sup>=</sup> <sup>∪</sup>*Nc <sup>i</sup>*=1*Ki*, where *K*<sup>i</sup> is coarse cell. The T*<sup>h</sup>* is a fine grid with *H* ≥ *h* ≥ 0, where *h* is size of fine grid. Each coarse-grid block (*Ki*) consists of the connected union of fine-grid blocks. Furthermore, a certain domain for this is constructed and denoted by *Nc* the coarse nodes number, by {*xi*}*Nc <sup>i</sup>*=<sup>1</sup> the vertices of the coarse mesh and define the neighborhood of the node *xi*:

$$
\omega\_i = \cup \{ \mathcal{K}\_i \in \mathcal{T}\_H \colon \mathfrak{x}\_i \in \overline{\mathcal{K}\_i} \}
$$

The following steps need to be implemented for GMsFEM realization:


In the first step of GMsFEM (Offline stage), the "snapshots" space must be constructed, a large dimensional snapshots space of local solutions. In the "snapshots" space, the following is considered:

$$-\operatorname{div}\left(k\_{\alpha}\nabla\psi\_{\vec{\jmath}}\right) = 0,\ \mathbf{x}\in\omega\_{i\nu}$$

$$\psi\_{\vec{\jmath}} = \delta\_{\vec{\jmath}}(\mathbf{x}),\ \mathbf{x}\in\partial\omega\_{i\nu}$$

where *k<sup>α</sup>* are the coefficients of the thermal conductivity, *δj*(*x*) are certain set of function defined on *∂ωi*, here *j* = 1, *Jω*. The *J<sup>ω</sup>* is a number of fine grid edges on *ω*. Therefore, following is defined:

$$V\_{\text{snap}} = \text{span}\{\boldsymbol{\psi}\_{j}^{\text{snap}} \colon 1 \le j \le J\_{\omega}\}, \text{ and } R\_{\text{snap}} = \left[\boldsymbol{\psi}\_{1}^{\text{snap}}, \dots, \boldsymbol{\psi}\_{J\_{\omega}}^{\text{snap}}\right].$$

This allow reducing the snapshot space to offline space via some spectral procedure. Offline space is constructed using the following local spectral problems in the snapshots space:

$$
\overline{A}\_{\text{off}} \,\overline{\Psi}\_k^{\text{off}} = \lambda\_k \overline{S}\_{\text{off}} \,\overline{\Psi}\_k^{\text{off}},
$$

where *A*off = *R*snapA*R<sup>T</sup>* snap, *S*off = *R*snapS*R<sup>T</sup>* snap, and

$$A = [a\_{\mathfrak{M}\mathfrak{n}}] = \int\_{\omega\_i^\circ} (k\_\mathfrak{a} \nabla \varrho\_{\mathfrak{m}\_\varGamma} \nabla \varrho\_\varPi) d\mathfrak{x}\_\varGamma \, S = [s\_{\mathfrak{m}\mathfrak{n}}] = \int\_{\omega\_i^\circ} (k\_\mathfrak{a} \varrho\_{\mathfrak{m}\_\varGamma} \varrho\_\varPi) d\mathfrak{x}\_\varGamma$$

To generate the offline space, the smallest *Mω<sup>i</sup>* off eigenvalues are chosen and the corresponding eigenvectors ψoff *<sup>k</sup>* = ∑ *m* Ψoff *mkϕ*off *<sup>m</sup>* are found for *<sup>k</sup>* = 1, 2, ... , *<sup>M</sup>ω<sup>i</sup>* off. The found eigenvectors must be multiplied by the partition of unity functions *χi*:

$$\varphi\_k^{\omega\_i} = \chi\_i \psi\_k^{\text{off}} \text{ for } 1 \le i \le N \text{ and } 1 \le k \le M\_{\text{off}'}^{\omega\_i}$$

where *Mω<sup>i</sup>* off denotes the number of eigenvectors that are sampled for each local *ωi*. After this procedure, conforming basis functions can be obtained in the space:

$$V\_{\rm off} = \text{span}\left\{\boldsymbol{\varrho}\_k^{\omega\_i}, \ 1 \le k \le M^{\omega\_i}, \ 1 \le i \le N\right\}.$$

Further, the projection matrix *R<sup>T</sup>* = 1 *ϕω<sup>i</sup>* <sup>1</sup> ,..., *<sup>ϕ</sup>ω<sup>i</sup> Mω<sup>i</sup>* 2 is defined. Using constructed multiscale space, the coarse-scale system is solved:

$$M\left(\frac{\partial T\_{\mathcal{L}}}{\partial t}\right) + A\_{\mathcal{L}}T\_{\mathcal{L}} = F\_{\mathcal{L}'} $$

where *Ac* = *RAf R<sup>T</sup>* and *Fc* = *RFf* .

After solving the coarse-scale solution, a solution on the fine grid *Tms* = *RTTc* can be calculated.

### **5. Numerical Results**

In this section, the effect of fibers in a pile made of different materials on the temperature field of the soil is calculated. Thus, the melting effect of the pile is numerically simulated.

Conditions for the impact of piles on the temperature regime of the soil:


The computational domain consists of several soil layers; it has one composite pile (Figure 2).

**Figure 2.** Computational domain (**a**) is structured, (**b**) is randomized.

There are 4 cases which were considered:


Soil temperature is 1.5 ◦C. Thermal-physical characteristics are taken from SP 25.13330.2012 and presented in Table 2. Calculations were performed for 365 days with a time step *τ* = 1 day (24 h).

**Table 2.** Thermal-physical characteristics of soils.


The temperature of the daylight surface of the structure base was set considering the amplitude of the air temperature fluctuations taken from the data of the Yakutsk meteorological station. Figure 3 demonstrates the changes in the air temperature. It shows the distribution of the air temperature during one year.

**Figure 3.** Air temperature.

Here are the results of numerical calculations. Figures 4 and 5 show the soil temperature distributions for different points in time. These results were obtained by the GMsFEM method for eight basis functions.

**Figure 4.** Temperature distribution for different time layers (*t* = 50, 200, 250, 365 days) when the fibers are structurally located. From left to right. (white line is isocline of zero of steel, black is isocline of zero of basalt).

**Figure 5.** Temperature distribution for different time layers (*t* = 50, 200, 250, 365 days) when the fibers are randomly located. From left to right. (white line is isocline of zero of steel, black is isocline of zero of basalt).

For numerical comparison of the fine-scale and multiscale solutions, relative *L*<sup>2</sup> and energy errors are used.

$$\|\|\left\|\left.\varepsilon\right\|\right\|\_{L\_2} = \sqrt{\frac{\int\_{\Omega} (T\_h - T\_{\rm ms})^2 d\mathbf{x}}{\int\_{\Omega} T\_h^2 d\mathbf{x}}},\\ \|\left\|\left.\varepsilon\right\|\right\|\_{a} = \sqrt{\frac{a\_{\Phi} \left(T\_h - T\_{\rm ms} T\_h - T\_{\rm ms}\right)}{a\_{\Phi} \left(T\_{\rm hr}, T\_{\rm hr}\right)}},\end{aligned}$$

where *Th* and *Tms* are the fine-scale and multiscale solutions.

Tables 3 and 4 show the relative errors of *L*<sup>2</sup> and energies for different number of multiscale basis functions and corresponding degrees of freedom (DOF). Degrees of freedom for fine scale mesh is denoted as DOF*f*. Here are the uncertainties for fibers arranged in structured order. When using basalt fibers, you can limit yourself to using 2 multiscale basis functions, and when using steel fibers, you need to use at least 4 multiscale basis functions.


**Table 3.** Relative *L*<sup>2</sup> and energy errors (%) for different number of multiscale basis functions (DOFf = 172,681). Case 1: structured basalt fibers.

**Table 4.** Relative *L*<sup>2</sup> and energy errors (%) for different number of multiscale basis functions (DOF*<sup>f</sup>* = 172,681). Case 2: structured steel fibers.


Tables 5 and 6 show the relative errors of *L*<sup>2</sup> and energies for different number of multiscale basis functions for case, when fibers are located randomly. It can be seen that using four multiscale basis functions leads to good accuracy.


**Table 5.** Relative *L*<sup>2</sup> and energy errors (%) for different number of multiscale basis functions (DOF*<sup>f</sup>* = 207,599). Case 3: random distribution of basalt fibers.

**Table 6.** Relative *L*<sup>2</sup> and energy errors (%) for different number of multiscale basis functions (DOF*<sup>f</sup>* = 207,599). Case 4: random distribution of steel fibers.


A comparison was also made for different fiber arrangements (Figures 6 and 7). As the figures show, random arrangement of the fibers is better. However, this may be due to the fact that they are located horizontally.

**Figure 6.** Temperature distribution for different time layers (*t* = 50, 200, 250, 365 days) in case of steel fiber reinforcement. From left to right. (red line is isocline of zero for structure locations, black is isocline of zero for random locations).

**Figure 7.** Temperature distribution for different time layers (*t* = 50, 200, 250, 365 days) in case of basalt fiber reinforcement. From left to right. (red line is isocline of zero for structure locations, black is isocline of zero for random locations).

### **6. Discussion**

For the solution of the Stefan problem related to the soil thawing effect of composite piles, GMsFEM was used. According to results showing the accuracy of the method, the usage of four basis functions leads to the results with error in L2 norm lower than 1%. The errors in energy norms are also very respectable. Thus, usage of GMsFEM for such problems is recommended.

Results show that the process of total freezing of piles with steel fiber inclusion in winter runs faster. Analysis shows that the length between isoclines of phase change is about 47 cm for structured and 37 for random distribution of fibers. Freezing of pile with steel fibers comes on the 239th day, and a pile with basalt inclusion freezes on the 250th day in the structured case. For random distribution, similar days are 245th and 252nd. Therefore, there are 11 and 7 days' delay between freezing of piles with steel and fiber inclusions, depending on the distribution.

### **7. Conclusions**

In this work, the effect of a composite pile on the temperature field of the surrounding frozen soil using GMsFEM is considered. Results of calculations using GMsFEM show good accuracy. It can be said that when using materials with higher thermal conductivity, the thermal regime restores faster. In the future, we plan to consider the problems in the 3D statement and work with the piles with more complex compositions.

**Author Contributions:** Conceptualization, P.V.S. and S.P.S.; methodology, P.V.S.; software, S.P.S.; investigation, P.V.S.; resources, P.S.; data curation, P.S.; writing—original draft preparation, P.V.S.; writing—review and editing, P.S.; visualization, S.P.S.; supervision, P.S.; project administration, P.V.S.; funding acquisition, P.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by grant RSCF 20-71-00133 and the mega-grant of the Russian Federation Government N14.Y26.31.0013.

**Data Availability Statement:** Data supporting reported results are available from the corresponding author by request.

**Acknowledgments:** The authors would like to thank every person/department who helped thorough out the research work. The careful review and constructive suggestions by the anonymous reviewers were gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

