In this research, we investigate at a heat and mass transfer process with artificial ground freezing. Almost every aspect of building and mining uses artificially frozen filter soils[andersland2003frozen, harris1995ground, newman2011artificial, alzoubi2017intermittent, jessberger1980theory]. For example, such technology is commonly utilized in mine sinking, tunneling, subway construction, and building construction, etc. Artificial freezing with cooling equipment is employed near piles in the construction of buildings on permafrost soils to assure stability by generating big chunks of frozen earth surrounding the pile, which will preserve the soil from defrosting throughout the summer [loveridge2012energy, fei2019ground, vasilyeva2017reduced]. This practice has been demonstrated to strengthen building foundations in the Far North.
The authors describe mathematical models of heat and mass transfer processes in frozen and thawed soils, as well as artificial freezing of filter soils in the works [esch2004thermal, alexiades2018mathematical]. A set of heat and mass transfer equations with a dynamic phase transition boundary is solved, which comprises a parabolic equation for temperature and an elliptic equation for pressure. We employ a finite element method to approximate problem on the fine grid. [bathe2007finite, szabo2021finite]. We apply the through counting approach to model temperature phase change [burago2021numerical, rubinshteuin1971stefan] and the fictional domain method to model pressure phase change [glowinski1994fictitious, glowinski1994fictitious1].
These processes are generally modeled in very huge domains where the construction body is located. It is important to construct a highly fine computational grid in order to precisely describe the process at each point in a large region. This method demands high computational prices and high-end computing hardware, such as a supercomputer. As a way out of the situation, You can utilize homogenization methods [talonov2016numerical, grigoriev2019effective, gavrilieva2018numerical, vasilyeva2017reduced] just design coarser grids to get out of the problem, but the accuracy of the solution will suffer. We propose a model reduction technique based on multiscale finite element method [efendiev2009multiscale].
Multiscale finite element method(MsFEM) is particularly suited for modeling problems in highly heterogeneous regions [hou1997multiscale, akkutlu2017multiscale]. Multiscale methods come in many forms, such as multiscale finite volume method (MsFVM)[hajibeygi2008iterative, lunati2006multiscale, lunati2008multiscale, tyrylgin2019embedded] which uses the finite volume method to generate the basis functions [eymard2000finite, moukalled2016finite]. Another MsFEM modification is the generalized finite element method(GMsFEM) [chung2016adaptive, chung2014generalized, efendiev2013generalized] which builds several bases in each local domain by solving local spectral problems. In the constraint energy minimizing generalized multiscale finite element method(CEM-GMsFEM), the basis building can be provided in an oversampled domain[chung2018constraint, cheung2018constraint, fu2020constraint]. A multiscale method that satisfies the properties of mass conservation is called mixed multiscale finite element method (Mixed MsFEM) [chen2003mixed, chung2015mixed, aarnes2008mixed]. An online generalized multiscale finite element method (Online GMsFEM) [akkutlu2016multiscale, chung2015residual, tyrylgin2021multiscale, spiridonov2021online] is particularly suited for nonlinear problems because it executes the procedure of enriching a multiscale space during the online stage of the method. A special type of multiscale basis functions based on constrained energy minimization problems are developed in [chung2018nonfrac, vasilyeva2019nonlocal, vasilyeva2019constrained, vasilyeva2019upscaling] and well-known as nonlocal multicontinuum method(NLMC).
We previously created a GMsFEM algorithm with an additional basis function for artificial ground freezing [vasilyeva2020multiscale]. We expand our technique and employ multiscale online basis functions to predict phase change in the heat and mass transfer problem with artificial ground freezing in this study. We perform a model reduction procedure which includes offline and online stage. We create offline basis functions using local spectral problems in the offline stage. To better forecast processes in mediums with high contrast, we use snapshot space in basis computation. The considered medium features layered heterogeneity with a large value jump between two layers. In the computational domain, our problem involves frozen pipes. We employ offline additional basis functions to account the effect of frozen pipes in order to make an accurate simulation. In the online stage, we construct an online multiscale basis functions. They assist in the accurate approximation of a shifting phase change boundary. Local residuals are used to compute online basis functions. Using online bases allows us to take smaller amount of offline basis functions with better accuracy.
The numerical results are shown in a two-dimensional heterogeneous domain. We consider a two type of boundary condition for pressure. In the first example we set flow from right to left and in the second example we set flow from top to bottom. We complete the validation process by showing the results for a variety of offline and online basis functions.
The paper organized as follows. In Section 2, we present a mathematical model for heat and mass transfer problem with numerical algorithm for phase change boundary. In Section 3, we present an approximation on the fine grid using finite element method. Next, in Section 4, we describe an algorithm of online generalized multiscale finite element method. Numerical results for two test cases are presented in Section 5.
2 Problem formulation
We consider the heat and mass transfer model with artificial ground freezing. Freezing pipes within the soils enable artificial ground freezing process. We use the classical Stefan model to simulate heat transfer processes with phase change [samarskiy2009computational, Rubinstein1967problem]. In this approximation, we assume that the phase change occurs at a given phase change temperature . Let is the domain of the liquid phase where the temperature exceeds the phase change temperature:
and is the domain of frozen phase:
We denote by , and , the density and specific heat capacity of the liquid and frozen zones, respectively. We define the indicator function as piecewise defined function
For the coefficients of heat capacity and thermal conductivity, we have
Thermal processes, which are accompanied by phase change, absorption and release of latent heat of fusion, are described by the equation
Here is the specific heat of the phase change, is the filtration speed in the soil.
We have heat capacity coefficients because we are considering heat transfer in a porous material:
here is the porosity and the subscripts , , and denote the skeleton of porous medium, water and ice, respectively. We use comparable terms for heat conductivity coefficients in the solid and frozen zones:
To take into account the filtration in the soil, we write the continuity equation for a completely saturated porous medium
We use Darcy’s law to express the relationship between the filtration rate and the pressure gradient:
is the absolute permeability tensor of the porous medium,is the fluid viscosity, is the free fall acceleration and is the density of water.
The system of equations (1), (4) is the base for modeling the processes of thermal stabilization of filter soils. In our problem the main difficulties of numerical simulation are generated by the phase change, we need to solve the filtration problem with a free (unknown) boundary. We use computational algorithms for through counting, which were used by many authors in solving similar problems of heat and mass transfer with phase transformations [samarskii1993numerical, belhamadia2012enhanced, bhattacharya2002fixed].
In our implementation we assume that the phase change occurs not when , but in some small interval near the temperature of phase change . Instead of the indicator function , we take its piecewise linear approximation of the following form:
Then, we have
Thus, instead of (1) in the computational domain , we solve the following equation for temperature:
This non-linear parabolic equation is supplemented with appropriate initial and boundary conditions.
Since our problem (4)-(5) has a moving phase transition boundary, we are faced with the problem of constructing a numerical algorithm for calculating the pressure. For the numerical solution of such a problem without rebuilding the computational grid, we use the method of fictitious domain, which is based on the transition to solving the problem in a wider domain. An approximate solution, depending on the continuation parameter , we will solve in the entire computational domain . When using a variant of the method of fictitious regions with continuation by leading coefficients, the solution is determined from the equation
Here, the discontinuous coefficient is determined by the expression
with very small . With this choice of coefficients, equation (6) simulates filtration in the region with a very small coefficient with the parameter .
As the boundary condition for the temperature, we take the zero Neumann boundary condition over the entire boundary of the computational domain:
In our problem we have an artificial ground freezing process. To simulate it, we set the following boundary condition on the pipes:
where is the boundary of artificial ground freezing pipes.
For the pressure, we set boundary conditions on the boundary of the entire computational domain :
When using the computational algorithm for through counting, it is necessary to pay special attention to setting the parameters and . These parameters indicate the width of the phase transition, they determine the accuracy of the phase change boundary for solving on a fine grid. and depend on the dimension of the grid, for an accurate approximation, a sufficient number of grid elements within the phase transition is required. In the paper, we chose these parameters based on the works [vabishevich2014numerical, vasil2020accurate] , where the authors studied the influence of hyperparameters on the solution of the problem with a phase change.
, where the authors studied the influence of hyperparameters on the solution of the problem with a phase change.
3 Fine grid approximation
Let us discretize in space using the finite element method of the initial-boundary value problem (5),(6), (8)-(11) for the resulting basic system of equations, which is used to simulate artificial freezing processes. To reformulate the original system of temperature and pressure equations, we multiply equation (5) by the test functions , and equation (6) by the test function . We integrate our system using Green’s formula and boundary conditions (6), (7) and obtain the following formulation:
Here is Sobolev space that consist of functions q such that and have a finite integral in and .
To approximate the equation for temperature in time, we use the standard implicit scheme. The simplest linearization is used when specifying the coefficients from the previous time layer.Let and be the time step. We denote the solution at time by . According to (12) we can derive a variational formulation:
Let write approximation in the matrix form:
Approximation in space (13) is associated with the choice of the corresponding finite element spaces for temperature and pressure and is carried out in the standard way. For both temperature and pressure, we use Lagrangian finite elements of the first order on a triangular mesh.
4 Coarse grid approximation
To make approximation on the coarse grid, we apply an Online Generalized Multiscale Finite Element method. A coarse grid is denoted by
here is the cell of coarse grid. Following that, we create a local coarse neighborhood domain it is produced by merging all coarse cells around one coarse grid vertex:
where and – the number of coarse grid vertices.
We can identify two stages in GMsFEM: online and offline:
Creating local coarse neighborhood domains from the coarse grid;
In each local domain we construct an offline multiscale basis functions;
Construction an additional basis function for temperature in each local domain that contains ground freezing pipes;
Assembling a projection matrix of the multiscale space. This mathix consist of offline basis functions.
Assembling a system of equations on fine grid and projecting it on the multiscale space by matrix;
Solving a system of equations on the coarse grid by the resulting approximation on a coarse grid;
Constructing an online multiscale basis functions on certain time layers;
Resolving system on the coarse grid using new matrix with online bases.
We should note, that we perform steps 3 and 4 only on certain time steps, otherwise we skip this steps. We add online bases not on every time layer, for example, we can do it on every 5-th or 10-th time layers. Following that, we go over the GMSFEM algorithm in further detail.
4.1 Offline stage
First, to build an offline multiscale basis functions, we need to define a snapshot space for temperature and for pressure . Therefore, we have to build projection matrices for these spaces. We solve the following problems in each local domain to identify the matrix’s members:
where is the discrete delta function which equal 1 at the -th fine grid node and zero elsewhere , is a number of fine grid vertices on boundary ). And is the outer boundary of local domain .
The purpose of creating a snapshot space must be discussed. Almost all applied problems are modeled in domains with high coefficient contrast. Such domains may contain fractures or channels. We will compute our problem in domain with high coefficient contrast. We also need to employ a snapshot space to account all important properties of the solution and produce a good approximation. The snapshot vectors can describe most important characteristics of computational domain and help to construct more accurate projection to multiscale space.
Now we can define a snapshot spaces and projection matrices for temperature and pressure:
The next step of offline multiscale basis construction is computing spectral problems in each local domain
. We should mention that we build a decoupled multiscale basis functions. To obtain a set of basis functions, we solve next eigenvalue problems:
Then we obtain solutions of spectral problems on space , . We generate offline multiscale basis functions using the initial eigenvectors that correspond to the smallest eigenvalues.
To get offline bases, we utilize the following formula:
Where is the linear partition of unity function that equals to standard coarse grid nodal basis function and . The example of first 4 offline multiscale basis functions in one local domain is presented in Figure 1.
After that, we’ll be able to create multiscale spaces
and projection matrices
Additional offline multiscale basis function.
In our problem, we have freezing pipes in modeling of temperature field. We must consider the impact of boundary condition . To calculate additional basis functions for temperature, we solve the next local problem in each local domain that includes freezing pipes
Then we multiply the solution of local problem (21) on partition of unity function to obtain an additional basis function . It should be noted, that the calculation of additional bases does not depend on the location of the freezing pipes, they can be both inside the coarse block and coincide with the coarse edges. One can see it in the solution of problem (21) in Figure 2.
As a result, we get the following multiscale space and projection matrix for temperature:
where is the number of with ground freezing pipes.
4.2 Online stage
One can see that in offline basis construction we didn’t take into account the phase change. We use thermal conductivity and porous medium permeability from liquid zone. We can’t take into account the phase change by offline bases, because it changes by time. We can neglect the phase transition, as we did earlier in [vasilyevastep2020multiscale], but in this paper we will take into account the phase change by online multiscale basis functions. The Online GMsFEM is well suited for solving non-linear problems. This method can significantly reduce the number of offline basis functions with good improvement in accuracy. It is better to add 1 or 2 online basis in each then computing a large number of offline bases. But this process has a big computational cost, because of this we compute online bases only in certain time layer, for example, in each 5-th or 10-th time layer. We construct an online multiscale basis function based on local residuals that provides fast error decay.
To make update of projection matrices in each -th time step, firstly, we need to solve the coarse grid system using only offline multiscale basis functions to find local residuals:
where , are multiscale solutions that reconstructed on fine grid, and are predefined projection matrices without enrichment with online bases at -th time layer
To complete online stage with solving coarse system with updating online multiscale bases functions, firstly, we solve the coarse scale system (22) using offline projection matrices and . Then we add online bases in these matrices and obtain new projection matrices , . We use these matrices until the next update procedure, we update projection matrices on each certain -th time layer. The process can be iterable to add more multiscale online basis functions:
where denotes the number of online iterations. On the time layer where we don’t update online bases, we use , with old online bases.
Next we present a construction of online multiscale basis functions that based on solution of the following problem in each local domain :
with the boundary condition on . After solving problem (23) we multiply the solution on partition of unity functions to obtain an online multiscale basis functions:
We demonstrate an example of online multiscale basis functions in Figure 3.
Next we update multiscale spaces with online basis functions:
Next, we use current solutions , to find a certain number of online basis functions in each local domain :
with , , .
5 Numerical results
In this paper, we consider problems in two-dimensional formulation in heterogeneous domain with high contrast. We consider modeling of heat and mass transfer in computational domain . For modeling, we take freezing pipes in computational domain. Our heterogeneity is represented by layers and simulates a layered reservoir. We construct structured triangular fine grid with 14641 vertices and 14400 cells. For multiscale solver, we generate coarse grid with 325 vertices and 288 cells. Computational domain, computational grid and schematic heterogeneity are presented in fig. 4.
We perform modeling of ground freezing process with days and times steps. In our simulation we set Dirichlet boundary condition on the freezing pipes and on the other boundaries we set . As initial condition for temperature we take . For pressure we take different boundary conditions:
Test 1. Flow from right to left, at the left boundary we set , on the right boundary – ,
Test 2. Flow from right to left, at the top boundary we set , on the bottom boundary – .
On the other boundaries for pressure, we apply .
We take the next values of coefficients:
We will examine the multiscale and fine grid solutions for numerical results validation. To do this, we use relative errors in and norms:
where is a solution on a fine grid obtained by finite element method, is a multiscale solution, instead of function we use temperature and pressure .
|Number of||Offline basis||1 Online basis||2 Online basis|
|Number of||Offline basis||1 Online basis||2 Online basis|