УДК 532.546: 551.340
Computer Simulation and Optimization of Soil Thawing Using Microwave Energy
Sergej A. Nekrasov* Vladimir S. VolkoV
Department of Applied Mathematics South-Russian State Polytechnical University Prosveschenia, 132, Novocherkassk, 346400
Russia
Received 19.10.2015, received in revised form 06.11.2015, accepted 15.12.2014 The article describes a mathematical model and numerical method for calculating and optimizing the speed of soil thawing using microwave energy. The main numerical method is the enthalpy method without smoothing concentrated heat. Examples of calculations and analysis of the results are presented. A series of calculations are carried out for the cases of sandstone and sandy loam thawing with the use of standard magnetrons. Nonlinear behavior of thermal parameters, phase transitions and boundary conditions significantly influence temperature distribution and electric field. The losses due to convective cooling and evaporation are negligibly small. The results can be used in research and in designing mining equipment for permafrost.
Keywords: thawing, soil, microwave, modeling, optimization. DOI: 10.17516/1997-1397-2016-9-1-60-68.
Introduction
The construction in permafrost zone is continuously increasing, both in Russia and in the World. In this regard, the most important problems in geological engineering and construction are forecast and control the temperature field distribution in soils with different lithological parameters [1-6]. Mathematical model and numerical procedure to forecast temperatures in a mine are presented in [7-9]. The model includes two-dimensional energy equation for air and Fourier equations to describe the thermal field in thawed and frozen rocks zones. These equations are complemented by Stephen equation on the moving phase boundary and Newton's equation on the boundary "air-rock". The initial conditions assume an arbitrary temperature distribution in space and time. Function that describes the temperature variation at the mine entrance takes into account not only seasonal but also daily fluctuations. It allows one to take into account climate features of the northern regions. This model also allows one to take into account the main factor affecting the intensity of heat transfer in time, namely, humidity (iciness) of rocks.
It was found [10-11] that changes in the properties of rocks during seasonal freezing are of great importance for the northern regions with harsh winters. Physical and mechanical properties of dispersed rocks in thawed and frozen state are studied for a long time. However, changes in the properties of rocks under the influence of freeze-thaw cycles as well as changes in the properties of rocks during freezing are poorly understood. Problems of the optimal control of temperature fields are relatively new and still little understood. One of these problems is to optimize the process of thawing of frozen soils using magnetrons. There are several approaches
* [email protected] [email protected] © Siberian Federal University. All rights reserved
to solve the problem. They include direct methods, variational methods and methods of optimal control. This article describes the use of a direct method of optimization. We consider the thawing technology that uses the microwave energy produced by the magnetrons. According to [1], the magnetron M-111 with operating frequency 915 MHz and excavating machine provide 0.3-0.5 m depth of excavating per pass. The energy consumption is about 5-6 kWh/m3. The magnetron M81 with operating frequency 2.45 GHz has an output of 5 kW. It provides the rate of penetration equal to 4-10 m/h with borehole diameter equal to 0.8 m in the case of frozen sandy loam. These parameters are approximate because energy consumption is proportional to temperature of frozen soil so it depends on the weather conditions.
Mathematical model
In the process of thawing there are the following stages of the thermal process: heating of frozen soil to a temperature of melting (solid phase), heating of unfrozen soil and further thawing of frozen soil (liquid phase), intense evaporation and even boiling of hot melt water can also take place.
Thermal processes in the soil are described in the context of non-stationary one-dimensional model. It is assumed that the soil occupies a half-space x> 0, x is the coordinate. Microwave energy is transmitted into the soil by means of transverse electromagnetic wave. The electric field strength on the soil surface (x=0) is equal to E0.
The thermal process is described by the multi-phase Stefan problem in the enthalpy formulation with concentrated heat capacity:
dH d2A r, ,
= dX2 + f (x,t),x> °,t> 0, (1)
d A
— = -q, x = 0, t> 0, (2)
dx
T = T0, x> 0,t = 0; T ^ T0, x , (3)
where H is the enthalpy
H (T) = / p (u) c(u)du + wf p (Tf i - 0) rf (T - Tf i) + wf p (Tf2 - 0) rf2$ (T - Tf2),
J To
t is time, $ (T) is the Heaviside function, T is the temperature, Tf i is the temperature of ice melting, rf i is the latent heat of melting, Tf2 is the boiling temperature of water, rf2 is the latent heat of vaporization of water, T0 is the initial soil temperature (T0<Tf i), p=p (T) is the density, pi (T) is the density of frozen soil (T<Tf i), p2 (T) is the density of the unfrozen soil (Tf i < T<Tf 2), c=c (T) is the heat capacity of the soil, wf is the relative soil moisture,
rT
A (T) = I A (u) du, where A (T) is the thermal conductivity, q is the density of the surface
To
source. Taking into account the convective heat loss and evaporation from the surface of the unfrozen soil, we have
q=k(T-Tair )+P2(T )rf 2v.
vap :
where k is the coefficient of convective heat transfer, rf2 is the specific heat of evaporation, vvap is the speed of the vaporization front. This speed is approximately determined by the
(Ts)
2nRTs p2 (Ts) '
corresponding to the surface temperature of the thawed soil, Ts = T(0,t), M is the molar mass
Hertz-Knudsen equation as vvap=\/ n_r>rTi n ^^, where pn is the saturated vapor pressure
of water and R is the gas constant. Assuming that heat of evaporation qvap is independent of temperature, the vapor pressure is derived from the Clausius-Clapeyron equation:
,m \ (MrVap MrVap \
pn (Ts) = -RT) '
I bp
where p0 is the external pressure, Tbp is the boiling temperature at pressure p0.
In equation (1) f is the specific heat release under microwave heating of dielectric:
//
f=U£r £0E2,
where w is the angular frequency of the electromagnetic field [12-15]. Complex permittivity is equal to
£r — £r i£r ,
where £r is the real part of complex permittivity (relative permittivity) and £r is the imaginary part of complex permittivity of the soil (dielectric loss factor), £0 is the dielectric constant, E is the electric field strength. Dielectric loss factor is £r — £ tgd. If the conductivity of a material a is small or the frequency is high so that a C w£ (£ = £r £0) then the dielectric heating is the dominant mechanism to transfer the energy of electromagnetic field (EMF) into material. The coefficient characterizing the depth of penetration of EMF into the soil is the imaginary part of the wave number: b = 3—^^ \J(\£r | — Re£r) /2.
Finite-difference method
Numerical method is based on the multi-phase Stefan problem in enthalpy formulation with concentrated heat capacity [2-6]. The problem domain is infinite so it is advisable to use a non-uniform grid. Equation (1) is approximated as follows
Hp+1—np_ 1 (Ap+1—Ap Ap —Ap_!
j—Ap — Atj±)
\ hx,j+1 I
r . h h 1 + fp, j = l,...,n — 1; p = 0,1,...,
hx,j y hx,j+1 hx,j I
where t is the time step value ( tp=pr), p is the number of time step, hx,j is the grid size,
v h^ x j \ h^ - 1 1
j is the number grid node, hx,j = —■——, Tjp is the temperature grid function, Tjp &
'2 j j T (xj ,tp) ,j=0,... ,n; p= 0,10...; Hp = H j ; Ap = A j.
The initial condition is Tj0 = T0, j = 0,... ,n. The first order approximation of boundary
condition with respect to h=maxj=1i...,n hx,j is
Tp Tp
a (T0p) Th= q (Tip) ,p =1, 2,....
hx,1
The second order approximation is
Tp= (hx,1 + hx,2 )2 tp__^xA_tp — hx,1 (hx,1 + hx,2) q (T0)
0 (2hx,1 + hx,2)hx,2 1 (2hx,1 + hx,2)hx,2 2 2hx,1 + hx,2 A (Tp) '
The Newton-Raphson method is used to solve the last two nonlinear equations. Condition on the right border of the grid is Hnp= 0,p = 0,1,....
The grid size is dictated by the requirement of accuracy of the finite-difference solution. Near the boundary the grid refinement is applied to improve the accuracy of approximation of the boundary condition.
The time integration scheme is explicit. Explicit schemes are preferred over implicit schemes. Explicit schemes are cost-effective and to a lesser degree they break the symmetry of the problem. This is of particular importance for nonlinear problems. The value of time step is chosen according to the well known stability condition of explicit schemes for the heat equation.
The employed enthalpy method does not use smoothing of concentrated specific heat. Difficulties emerge when concentrated heat capacity is taken into account [4].
Similar problems of modeling of phase transitions in soils are discussed in [16].
Soil properties
Dry soil has high resistivity and relatively low dielectric constant, for example, the average sand resistivity is equal to 3000 Ohm*m and the relative dielectric constant is around 4. The presence of water in the soil greatly increases its conductivity. Wet sand has resistivity of about 100 Om*m. The soil often contains a large amount of bound water and in the frozen soil a portion of water remains unfrozen. This also reduces the resistance to 30 Om*m [12-15]. Thus, the electrical properties of water are important. Fig. 1 and 2 show the relative dielectric constant of water and ice. Dielectric loss in the soil is mainly due to dipole polarization and electric conductivity. The dielectric constants for the widely used frequency 915 MHz are: for the frozen soil e'r = 4(1 — wf) + 3wf, e'r = 0.0004(1 — wf) + 0.0003wf; for the unfrozen soil e'r = 4(1 — Wf) + 70wf; e'r = 0.0004(1 — Wf) + 55wf, where Wf is the soil moisture.
Fig. 1 and 2 show that often used magnetron frequencies 0.5-2 GHz are best suited for heating water and unfrozen soil (in both cases the depth of penetration of EMF is about 1 cm) but they are relatively less effective for heating ice and frozen soil (penetration depth of EMF is about 0.5 m). Convective heat transfer coefficient k is taken to be 1W/(m2*K), soil moisture wf =0,1. Air temperature and the initial temperature of the soil are assumed to be equal. They vary between -50 and 0 oC. Thermophysical parameters of various types of soil as a function of moisture value are presented in [17]. Typical values for 12.
Fig. 1. The dependence of the real part of the relative permittivity of ice water and the frequency of the electromagnetic waves
Fig. 2. The dependence of the imaginary part of the relative permittivity of ice water and the frequency of the electromagnetic waves
Magnetron characteristics
Normal power of magnetrons used in construction works are P = 5 kW and efficiency is about 70%. Radiated specific instantaneous power of magnetron is the magnitude of the Poynting vector P = EH and the actual value <P> = 0.5 EH, where E is the electric field strength and H is the magnetic field strength. Since H = E/Z, where Z is the impedance of vacuum (about 377 ohms), the average radiated power is < Prad >= 0.5E2S/Z, where S is the surface area of radiation. Assuming that efficiency is 70% (< Prad >= 0.7P) and the surface area of radiation is S = 0.5m2, we find the following estimate of the electric field strength emitted by the magnetron: E = (1.4PZ/S)0 5 = (1.4 * 5000 * 377/0.5)0 5 « 2300 V/m. There are often higher flux densities of microwave energy with electric field strength is between 7000 and 11 000 V/m [1].
Simulation results
Data used in simulations: the size of the region Lx =5 m, the number of grid nodes n = 50 — 2000; EMF frequency f =915 MHz; conductivity of phases (sandy loam soil): a1= 1/150 S/m, ^2= 1/25 S/m, as = 1/150 S/m; electric field strength Eo= 11000 V/m; borehole radius R0 = 0.4 m; the required depth of thawing h = 0.25 m. The characteristic depth (in meters) of penetration of microwave energy for each of the phases (sandy loam soil) for 915 MHz: Li = 0.65; L2 = 0.03; L3 = 0.80; for 430 MHz: Li = 0.72; L2 = 0.05; L3 = 0.80. After the start of melting the flux of microwave energy in the solid phase is reduced because the field is strongly attenuated in the liquid phase.
Fig. 3-5 show the space distributions of the temperature field (time interval is 10 sec). The duration of exposure to microwave radiation is about 50 sec. Fig. 3-5 show that temperature distribution depends on the initial temperature. In the case shown in Fig. 3 thawing time is about 594 sec, volumetric source energy is 32 kWh or 251 kWh per cubic meter. In the case shown in Fig. 4 thawing time is about 422 sec, volumetric source energy is 24 kWh or 190 kWh
per cubic meter. In both cases one needs 2.1 kWh (16.8 kWh per cubic meter) to melt ice only.
Characteristics of the thawing process (energy consumption, time) also significantly depend on the soil moisture. In the case shown in Fig. 5 thawing time for the soil moisture 2 % is approximately 380 sec, volumetric source energy is 18 kWh or 142 kWh per cubic meter. In both cases one needs 0.4 kWh (3.25 kWh per cubic meter) to melt ice only, xmeit(tF) = 25 cm, xbp(tF)= 16 cm, the surface temperature is 99°C and the maximum temperature is 108 °C.
Thermo-physical processes in the case of high-power microwave heating of moist soil are similar to the processes in microwave ovens in the case of food defrost. There are phase transitions ice-water and water-vapour.
Fig. 3. Space distribution of the temperature field T0 = — 30°C, wf = 0.1
Fig. 4. Space distribution of the temperature field T0 = —5°C, wf = 0.1
Error of numerical solution is estimated by the heat balance (it is about 1%). Fig. 6 shows positions (in arbitrary units) of phase transitions fronts versus time (To = —5 °C,Wf = 0,1). The maximum values of the coordinates of the phase transition fronts are
xmeit(tF)=25 cm and xbp(tF)=20 cm. The surface temperature is 99 °C and the maximum temperature is 100 °C.
Fig. 5. Space distribution of the temperature field T0 = —30 °C, wf = 0.02
Fig. 6. Line 1 is temperature of the surface, lines 2 and 3 are coordinates of front of melting and front of boiling T0 = -5 °C
Optimization
Due to the errors introduced by the finite-difference method, the objective function can have the multitude of false and weak extrema. It is advisable to carry out preliminary smoothing either to use statistical methods of optimization or direct search method. The objective function with respect to the thawing time to a predetermined depth and energy consumption is monotonic: energy consumption increases and the thawing time decreases with increasing of the electric field strength E0. Consider the objective function with respect to the depth of thawing h and penetration H for a working day (8 hours). The objective function is g(h)=8*3600*h/t(h). The objective function is shown in Fig. 7. Optimal values are hopt= 0,39 m and Hopt= 17,86 m.
g, m/day
20 Hop! 16
12
8
4
----—i—-----—»
0 16 32 liopi 48 64 80li, cm
Fig. 7. Objective function g(h)
Conclusions
A mathematical model and numerical method for calculating and optimizing the speed of soil thawing using magnetrons are presented. The enthalpy method without smoothing concentrated heat is used for simulations. A series of calculations are carried out for the cases of sandstone and sandy loam thawing with the use of standard magnetrons. It is found that nonlinear behavior of thermal parameters, phase transitions and boundary conditions significantly influence temperature distribution and electric field. It is also found that the losses due to convective cooling and evaporation are negligibly small. The results can be used in research and in designing mining equipment for permafrost.
References
[1] V.M.Petrov, Softening rocks powerful electromagnetic field of the microwave, Inf. Ra-dioelektronika i Kommunikatsii, 22(2002), no. 4, 63-73 (http://www.rit.informost.ru/rit/4-2002/63.pdf) (in Russian).
[2] S.Melnikov, Mathematical modeling of time-dependent temperature field in two-phase media, Nauka i obrazovanie, no. 2, 2012 (http://technomag.edu.ru/doc/330390.html) (in Rus-
[3] I.A.Gishkelyuk, Yu.V.Stanilovskaya, Computer 3D modeling halo thawing soils with ice wedges around the pipeline, Truboprovodny transport: teoriya i praktika, 40(2013), no. 6, 14-20 (in Russian).
[4] N.A.Buchko, Enthalpy method of numerical solution of heat conduction problems in freezing and thawing soils, SPbGUNTiPT, 2015 (http://refportal.com/upload/files /entalpi-iny_metod_chislennogo_resheniya.pdf) (in Russian).
[5] A.R.Pavlov, M.V.Matveeva. Iterative finite-difference scheme for the problem of heat and mass transfer during the freezing of soils, Vestnik of SSU - natural science series, 56(2007), no. 6, 242-253 (http://vestnik-samgu.samsu.ru/est/2007web6/math/2007560310.pdf) (in Russian).
sian).
[6] B.G.Aksenov et al., Simulation of oscillations of the boundary of freezing and thawing in the ground and external building structures, Sovremennye problemy nauki i obrazovaniya, (2012), no. 2 (URL: www.science-education.ru/102-6088) (in Russian).
[7] A.F.Galkin, Unsteady heat transfer coefficient for the bottom zone of mining mines of the North, IPTPN SB RAS, Yakutsk, 1978, 100-105 (in Russian).
[8] A.F.Galkin, I.N.Los, Assessing the impact of geothermic deposits in the choice of strategy working out mine field, Translated from Fiziko-Tekhnicheskie Problemy Razrabotki Poleznykh Iskopaemyh (FTPRPI), Novosibirsk, 1985, no. 2, 86-89 (in Russian).
[9] A.F.Galkin, Prognosis and selection of the optimal parameters of the thermal regime in the construction, operation and integrated use of mines in permafrost. Abstract of dissertation for the degree of Doctor of Technical Sciences, Saint Petersburg, 2009.
[10] A.S.Kurilko, Control of physical and mechanical properties of rocks at alternating temperature influence, Abstract of dissertation for the degree of Doctor of Technical Sciences, Ekaterinburg, 2005 (in Russian).
[11] A.S.Kurilko, M.D.Novopashin, Investigation of the effect of temperature on the strength of rocks, Translated from Fiziko-Tekhnicheskie Problemy Razrabotki Poleznykh Iskopaemyh (FTPRPI), 2005, no. 2, 32-36.
[12] The electrical conductivity and polarizability of permafrost in a constant electric field. Library materials engineering, 2014 (https://injzashita.com/elektroprovodnost-i-polyarizuemost-merzlix-porod-v-postoyannom-elektricheskom-pole.html) (in Russian).
[13] The values of the calculated resistivity of soil. Table, 2015 (http://www.zandz.ru/udelnoe_ soprotivlenie_grunta.html) (in Russian).
[14] A.N.Shuvaev, D.A.Genze, The dielectric constant of the soil broken structure, Vestnik TGASU, (2011), no. 1, 200-206 (http://www.tsuab.ru/upload/files/ additional/ 20_SHu-vaev_file_3891_3491_2257.pdf) (in Russian).
[15] The dielectric properties of water and ice, Laboratory meteotehnology, 2015 (http:// www.meteolab.ru/projects/dielectric/) (in Russian).
[16] V.I.Vasilev, A.M.Maksimov, E.E.Petrov, G.G.Tsypkin, Mathematical model of the freezing-thawing of saline frozen soil, Journal of Applied Mechanics and Technical Physics, 36(1995), no. 5, 689-696.
[17] Thermal properties of soils (electronic resource), Electron. dan. - Reference and information portal, 2012, Access: http://minkor.ru/upload/spravochnik /190609-2.pdf, free. - Caps. from the screen.
Компьютерное моделирование и оптимизация оттаивания грунтов при помощи энергии СВЧ
Сергей А. Некрасов Владимир С.Волков
В статье рассмотрены математическая модель и численный метод для расчета и оптимизации
процесса скоростного оттаивания грунтов при помощи энергии СВЧ. Приводятся соответствующие примеры расчетов и анализ результатов.
Ключевые слова: оттаивание, грунт, СВЧ, моделирование, оптимизация.