A convection-conduction model for analysis of the freeze-thawconditions in the surrounding rock wall of atunnel in permafrost regionsHE Chunxiong( 何春雄 )，(State Key Laboratory of Frozen Soil Engineering, Lanzhou Institute of Glaciology andGeocryology,Chinese Academy of Sciences, Lanzhou 730000, China; Department of Applied Mathematics, South China University of Technology, Guangzhou 510640, China)WU Ziwang( 吴紫汪 )and ZHU Linnan( 朱林楠 )(State key Laboratory of Frozen Soil Engineering, Lanzhou Institute of Glaciology andGeocryologyChinese Academy of Sciences, Lanzhou 730000, China)Received February 8, 1999AbstractBased on the analyses of fundamental meteorological and hydrogeological conditions at the site of a tunnel in the cold regions, a combinedconvection-conduction model for air flow in the tunnel and temperature field in the surrounding has been constructed. Using the model, the air temperature distribution in the Xiluoqi No. 2 Tunnel has been simulated numerically. The simulated results are in agreement with the data observed. Then, based on the in situ conditions of sir temperature, atmospheric pressure, wind force, hydrogeology and engineering geology, the air-temperature relationship between the temperature on the surface of the tunnel wall and the air temperature at the entry and exit of the tunnel has been obtained, and the freeze-thaw conditions at the Dabanshan Tunnel which is now under construction is predicted.Keywords: tunnel in cold regions, convective heat exchange and conduction, freeze-thaw.A number of highway and railway tunnels have been constructed in the permafrost regions and their neighboring areas in China. Since the hydrological and thermalconditions changed after a tunnel was excavated， the surrounding wall rock materials often froze, the frost heaving caused damage to the liner layers and seeping water froze into ice diamonds，which seriously interfered with the communication and transportation. Similar problems of the freezing damage in the tunnels also appeared in other countries like Russia, Norway and Japan .Hence it is urgent to predict the freeze-thaw conditions in the surrounding rock materials and provide a basis for the design，construction and maintenance of new tunnels in cold regions.Many tunnels，constructed in cold regions or their neighbouring areas，passthrough the part beneath the permafrost base .After a tunnel is excavated，the original thermodynamical conditions in the surroundings are and thaw destroyed and replacedmainly by the air connections without the heat radiation, the conditions determined principally by the temperature and velocity of air flow in the tunnel ， the coefficients of convective heat transfer on the tunnel wall， and the geothermal heat. In order to analyze and predict the freeze and thaw conditions of the surrounding wall rock of a tunnel，presuming the axial variations of air flow temperature and the coefficients of convective heat transfer, Lunardini discussed the freeze and thaw conditions by the approximate formulae obtained by Sham-sundar in study of freezing outside a circular tube with axial variations of coolant temperature .We simulated the temperature conditions on the surface of a tunnel wall varying similarly to the periodic changes of the outside air temperature .In fact，the temperatures of the air and the surrounding wall rock material affect each other so we cannot find the temperature variations of the air flow in advance; furthermore，it is difficult to quantify the coefficient of convective heat exchange at the surface of the tunnel wall .Therefore it is not practicable to define the temperature on the surface of the tunnel wall according to the outside air temperature .In this paper, we combine the air flow convective heat ex-change and heat conduction in the surrounding rock material into one model， and simulate the freeze-thaw conditions of the surrounding rock material based on the in situ conditions of air temperature，atmospheric pressure，wind force at the entry and exit of the tunnel， and the conditions of hydrogeology and engineering geology. Mathematical modelIn order to construct an appropriate model, we need the in situ fundamental conditions as a ba-sis .Here we use the conditions at the scene of the Dabanshan Tunnel. The Dabanshan Tunnel is lo-toted on the highway from Xining to Zhangye, south of the Datong River,at an elevation of 3754.78-3 801.23 m, with a length of 1 530 m and an alignment from southwest to northeast. The tunnel runs from the southwest to the northeast.Since the monthly-average air temperature is beneath 0`}C for eight months at the tunnel site each year and the construction would last for several years，the surrounding rock materials would become cooler during the construction .Weconclude that, after excavation, the pattern of air flow would depend mainly on the dominant wind speed at the entry and exit ， and the effects of the temperature difference between the inside and outside of the tunnel would be very small .Since the dominant wind direction is northeast at the tunnel site in winter, the air flow in the tunnel would go from the exit to the entry. Even though the dominant wind trend is southeastly in summer, considering the pressure difference, the temperature difference and the topography of the entry and exit，the air flow in the tunnel would also be from the exit to entry .Additionally ，since the wind speed at the tunnel site is low，we could consider that the air flow would be principally laminar.Based on the reasons mentioned，we simplify the tunnel to a round tube，and consider that theair flow and temperature are symmetrical about the axis of the tunnel ，Ignoring the influence of the air temperature on the speed of air flow, we obtain the following equation:where t，x，r are the time，axial and radial coordinates; U，V are axial and radial wind speeds; T is temperature; p is the effective pressure(that，isair pressure divided by air density); v is the kinematic viscosity of air; a is the thermal conductivity of air; L is the length of the tunnel; R is the equivalent radius of the tunnel section; D is the length of time after the tunnel construction;，S f(t),S u(t)are frozen and thawed parts in the surrounding rock materials respectively; f ,u and C f, C u are thermal conductivities and volumetric thermalcapacities in frozen and thawed parts respectively; X= (x , r) ， (t) is phase change front; Lh is heat latent of freezing water; and To is critical freezing temperature of rock ( here we assume To= -0.1℃).2 used for solving the modelEquation(1)shows flow. We first solve those concerning temperature at that the temperature of the surrounding rock does not affect the speed of air equations concerning the speed of air flow, and then solve those equations every time elapse.2. 1 Procedure used for solving the continuity and momentum equationsSince the first three equations in(1) are not independent we derive the secondequation by xand the third equation by r. After preliminary calculation we obtain the following elliptic equation concerning the effective pressure p:Then we solve equations in(1) using the following procedures:(i ) Assume the values for U0， V0;( ii ) substituting U0 ，V0 into eq. (2)，and solving (2)，we obtain p0;(iii)solving the first and second equations of(1)，we obtain U0， V1;(iv)solving the first and third equations of(1)，we obtain U2，V2;(v)calculating the momentum-average of U1，v1 and U2，v2，we obtain the new U0，V0;then return to (ii);(vi)iterating as above until the disparity of those solutions in two consecutiveiterations is sufficiently small or is satisfied，we then take those values of p0，U0 and V0 as the initial values for the next elapse and solve those equations concerning the temperature..2 .2Entire method used for solving the energy equationsAs mentioned previously， the temperature field of the surrounding rock and the air flow affect each other. Thus the surface of the tunnel wall is both the boundary of the temperature field in the surrounding rock and the boundary of the temperature field in air flow .Therefore ， it is difficult to separately identify the temperature on the tunnel wall surface ，and we cannot independently solve those equations concerning the temperature of air flow and those equations concerning the temperature of the surrounding rock .In order to cope with this problem，we simultaneously solve the two groups of equations based on the fact that at the tunnel wall surface both temperatures are equal .We should bear in mind the phase change while solving those equations concerning the temperature of the surrounding rock，and the convection while solving those equations concerning the temperature of the air flow, and we only need to smooth those relative parameters at the tunnel wall surface .The solving methods forthe equations with the phase change are the same as in reference [3].2.3Determination of thermal parameters and initial and boundary conditions2.3.1 Determination of the thermal parameters. Using p= 1013.25-0.1088 H ，we calculatePGT where T is the yearly-average absolute air temperature， and G is the humidity constant of air. Letting C be the thermal capacity with fixed pressure,thePthermal conductivity，and the dynamic viscosity of air flow, we calculate the thermal conductivity andkinematic viscosity using the formulas ａand. The thermal parametersC Pof the surrounding rock are determined from the tunnel site.2 .3.2 Determination of the initial and boundary conditions .Choose the observed monthly average wind speed at the entry and exit as boundary conditions of wind speed， and choose the relative effective pressure p=0 at the exit ( that，isthe entry of the dominant wind trend) and p (1 kL / d ) v2/ 2[ 5] on the section of entry ( thatis，the exit of the dominant wind trend )，where k is the coefficient of resistance along the tunnel wall, d = 2R ，and v is the axial average speed. We approximate T varying by the sine law according to the data observed at the scene and provide a suitable boundary value based on the position of the permafrost base and the geothermal gradient of the thaw rock materials beneath the permafrost base.3 A simulated exampleUsing the model and the solving method mentioned above，we simulate the varying law of the air temperature in the tunnel along with the temperature at the entry and exit of the Xiluoqi No.2 Tunnel .We observe that the simulated results are close to the data observed[6].The Xiluoqi No .2 Tunnel is located on the Nongling railway in northeasternChina and passes through the part beneath the permafrost base .It has a length of 1160 m running from the northwest to the southeast, with the entry of the tunnel in the northwest， and the elevation is about 700 m. The dominant wind direction in the tunnel is from northwest to southeast, with a maximum monthly-average speed of 3m/s and a minimum monthly-average speed of 1 .7 m/s . Based on the data observed，we approximate the varying sine law of air temperature at the entry and exit with yearly averages of -5℃， -6.4℃ and amplitudes of 18.9℃ and 17.6℃ respectively. The equivalent diameter is 5 .8m，and the resistant coefficient along the tunnel wallis 0.025.Since the effect of the thermal parameter of the surrounding rock on the air flow is much smaller than that of wind speed ，pressure and temperature at the entry and exit ， we refer to the data observed in the Dabanshan Tunnel for the thermal parameters.Figure 1 shows the simulated yearly-average air temperature inside and at the entry and exit of the tunnel compared with the data observed .We observe that the difference is less than 0 .2 `C from the entry to exit.Figure 2 shows a comparison of the simulated and observed monthly-averageair temperature in-side (distance greater than 100 m from the entry and exit) the tunnel. We observe that the principal law is almost the same， and the main reasonfor the difference is the errors that came from approximating the varying sine law atthe entry and exit; especially , the maximum monthly-average air temperature of1979 was not for July but for August.4Prediction of the freeze-thaw conditions for the Dabanshan Tunnel4.1 Thermal parameter and initial and boundary conditionsUsing the elevation of 3 800 m and the yearly-average air temperature of -3℃, wecalculate the air density p=0 .774 kg/m 3 .Since steam exists In the air, we choose thethermal capacity with a fixed pressure of air C p 1.8744kJ /( kg.0 C), heatconductivity 2.0 10 2 W /( m.0 C ) andand the dynamic viscosity9.218 10 6kg /(m.s). After calculation we obtain thethermal diffusivity a= 1 .3788 10 5m 2 / s and the kinematic viscosity ，1.19 10 5 m 2 / s .Considering that the section of automobiles is much smaller than that of thetunnel and the auto-mobiles pass through the tunnel at a low speed ，we ignore the piston effects ，coming from the movement of automobiles ，in the diffusion of the air.We consider the rock as a whole component and choose the dry volumetric cavity d2400kg / m 3 ,content of water and unfrozen water W=3% and W=1%, and thethermal conductivityu1.9W / m.o c , f2.0W / m.o c ,heat capacityC V 0.8kJ / kg.o c and C f(0.8 4.128w u ) d , C u (0.8 4.128w u ) d1 W 1 WAccording to the data observed at the tunnel site ，the maximum monthly-average wind speed is about 3 .5 m/s ， and the minimum monthly-average wind speed is about 2 .5 m/s .We approximate the wind speed at the entry and exit asv(t ) [0.028 (t 7)2 2.5]( m / s) , where t is in month. The initial wind speed in thetunnel is set to beU ( 0, x , r ) U a (1( r ) 2 ), V ( 0 , x , r ) 0 . RThe initial and boundary values of temperature T are set to bewhere f(x) is the distance from the vault to the permafrost base ，and R0=25 m is the radius of do-main of solution T. We assume that the geothermal gradient is 3%，the yearly-average air temperature outside tunnel the is A=-3 0C ， and the amplitude isB=12 0C .As for the boundary of R=Ro，we first solve the equations considering R=Ro as the first type of boundary; that is we assume that T=f(x) 3% 0C on R=Ro. We find that, after one year, the heat flow trend will have changed in the range of radius between 5 and 25m in the surrounding rock.. Considering that the rock will be cooler hereafter and it will be affected yet by geothermal heat, we appoximately assume that the boundary R=Ro is the second type of boundary; that is，we assume that the gradient value，obtained from the calculation up to the end of the first year after excavation under the first type of boundary value, is the gradient on R=Ro of T.Considering the surrounding rock to be cooler during the period of construction，we calculatefrom January and iterate some elapses of time under the same boundary. Then welet the boundaryvalues vary and solve the equations step by step(it can be proved that the solution will not depend on the choice of initial values after many time elapses ).4 .2Calculated resultsFigures 3 and 4 show the variations of the monthly-average temperatures on the surface of the tunnel wall along with the variations at the entry and exit .Figs .5 and 6 show the year when permafrost begins to form and the maximum thawed depth after permafrost formed in different surrounding sections.4 .3Preliminary conclusionBased on the initial-boundary conditions and thermal parametersmentioned above, we obtain the following preliminary conclusions:1)The yearly-average temperature on the surface wall of the tunnel is approximately equal to the air temperature at the entry and exit. It is warmer duringthe cold season and cooler during the warm season in the internal part (more than100 m from the entry and exit) of the tunnel than at the entry and exit . Fig .1 showsthat the internal monthly-average temperature on the surface of the tunnel wall is 1.2℃higher in January, February and December, 1℃ higher in March and October, and1 .6℃ lower in June and August, and 2qC lower in July than the air temperature at the entry and exit. In other months the infernal temperature on the surface of the tunnel wall approximately equals the air temperature at the entry and exit.2)Since it is affected by the geothermal heat in the internal surrounding section，especially in the central part, the internal amplitude of the yearly-averagetemperature on the surface of the tunnel wall decreases and is 1 .℃6 lower than that at the entry and exit.3 ) Under the conditions that the surrounding rock is compact , without a great amount of under-ground water, and using a thermal insulating layer(as designed PU with depth of 0.05 m and heat conductivity=0.0216 W/m℃， FBT with depth of 0.085 m and heat conductivity =0.0517W/m℃ )，in the third year after tunnel construction，the surrounding rock will begin to form permafrost in the range of 200 m from the entry and exit .In the first and the second year after construction, the surrounding rock will begin to form permafrost in the range of 40 and 100m from the entry and exit respectively .In the central part ， more than 200m from the entry and exit, permafrost will begin to form in the eighth year. Near the center of the tunnel ，permafrost will appear in the 14-15th years. During the first and second years after permafrost formed，the maximum of annual thawed depth is large (especially in the central part of the surrounding rock section) and thereafter it decreasesevery year. The maximum of annual thawed depth will be stable until the 19-20th years and will remain in s range of 2-3 m.4)If permafrost forms entirely in the surrounding rock ，the permafrost will provide a water-isolating layer and be favourable for communication and transportation .However, in the process of construction ， we found a lot of underground water in some sections of the surrounding rock .It will permanently exist in those sections，seeping out water and resulting in freezing damage to the liner layer. Further work will be reported elsewhere.严寒地区隧道围岩冻融状况分析的导热与对流换热模型何春雄吴紫汪朱林楠（中国科学院寒区旱区环境与工程研究所冻土工程国家重点实验室）（华南理工大学应用数学系）摘要通过对严寒地区隧道现场基本气象条件的分析，建立了隧道内空气与围岩对流换热及固体导热的综合模型;用此模型对大兴安岭西罗奇 2 号隧道的洞内气温分布进行了模拟计算，结果与实测值基本一致;分析预报了正在开凿的祁连山区大坂山隧道开通运营后洞内温度及围岩冻结、融化状况.关键词严寒地区隧道导热与对流换热冻结与融化在我国多年冻土分布及邻近地区，修筑了公路和铁路隧道几十座.由于隧道开通后洞内水热条件的变化;，普遍引起洞内围岩冻结，造成对衬砌层的冻胀破坏以及洞内渗水冻结成冰凌等，严重影响了正常交通 .类似隧道冻害问题同样出现在其他国家（苏联、挪威、日本等 )的寒冷地区 .如何预测分析隧道开挖后围岩的冻结状况，为严寒地区隧道建设的设计、施工及维护提供依据，这是一个亟待解决的重要课题 .在多年冻土及其临近地区修筑的隧道，多数除进出口部分外从多年冻土下限以下岩层穿过 .隧道贯通后，围岩内原有的稳定热力学条件遭到破坏，代之以阻断热辐射、开放通风对流为特征的新的热力系统 .隧道开通运营后，围岩的冻融特性将主要由流经洞内的气流的温度、速度、气—固交界面的换热以及地热梯度所确定 .为分析预测隧道开通后围岩的冻融特性， Lu-nardini 借用 Shamsundar研究圆形制冷管周围土体冻融特性时所得的近似公式，讨论过围岩的冻融特性.我们也曾就壁面温度随气温周期性变化的情况，分析计算了隧道围岩的温度场[3]. 但实际情况下，围岩与气体的温度场相互作用，隧道内气体温度的变化规律无法预先知道，加之洞壁表面的换热系数在技术上很难测定，从而由气温的变化确定壁面温度的变化难以实现 .本文通过气一固祸合的办法，把气体、固体的换热和导热作为整体来处理，从洞口气温、风速和空气湿度、压力及围岩的水热物理参数等基本数据出发，计算出围岩的温度场 .1数学模型为确定合适的数学模型，须以现场的基本情况为依据 .这里我们以青海祁连山区大坂山公路隧道的基本情况为背景来加以说明 .大坂山隧道位于西宁一张业公路大河以南，海拔 3754.78~3801.23 m，全长 1530 m ,隧道近西南—东北走向 .由于大坂山地区隧道施工现场平均气温为负温的时间每年约长8 个月，加之施工时间持续数年，围岩在施土过程中己经预冷，所以隧道开通运营后，洞内气体流动的形态主要由进出口的主导风速所确定，而受洞内围岩地温与洞外气温的温度压差的影响较小；冬季祁连山区盛行西北风，气流将从隧道出曰流向进口端，夏季虽然祁连山区盛行东偏南风，但考虑到洞口两端气压差、温度压差以及进出口地形等因素，洞内气流仍将由出口北端流向进口端 .另外，由于现场年平均风速不大，可以认为洞内气体将以层流为主基于以上基本情况，我们将隧道简化成圆筒，并认为气流、温度等关十隧道中心线轴对称，忽略气体温度的变化对其流速的影响，可有如下的方程:其中 t 为时间， x 为轴向坐标， r 为径向坐标 ;U, V 分别为轴向和径向速度， T 为温度， P 为有效压力 (即空气压力与空气密度之比少， V 为空气运动粘性系数， a 为空气的导温系数， L 为隧道长度， R 为隧道的当量半径， D 为时间长度S f(t ) ,S u (t ) 分别为围岩的冻、融区域.f, u分别为冻、融状态下的热传导系数，C f , C u 分别为冻、融状态下的体积热容量，X=(x,r) ,(t ) 为冻、融相变界面,To 为岩石冻结临界温度( 这里具体计算时取To=-0.100C )，L h为水的相变潜热 .2求解过程由方程 (1)知，围岩的温度的高低不影响气体的流动速度，所以我们可先解出速度，再解温度 .2.1 连续性方程和动量方程的求解由于方程 ((1)的前 3 个方程不是相互独立的，通过将动量方程分别对求导，经整理化简，我们得到关于压力P 的如下椭圆型方程 :x 和r于是，对方程 (1)中的连续性方程和动量方程的求解，我们按如下步骤进行:(1)设定速度U0 ,V0 ;( 2)将U0 ,V0代入方程并求解，得 P 0(3)联立方程 (1)的第一个和第二个方程，解得一组解U 1,V1;(4)联立方程 ((1)的第一个和第三个方程，解得一组解U2,V2;(5)对 ((3) ,(4)得到的速度进行动量平均,得新的U0 ,V0返回 (2) ;(6)按上述方法进行迭代，直到前后两次的速度值之差足够小.以 P0 , U0 , V0作为本时段的解 ,下一时段求解时以此作为迭代初值.2.2 能量方程的整体解法如前所述，围岩与空气的温度场相互作用，壁面既是气体温度场的边界，又是固体温度场的边界，壁面的温度值难以确定，我们无法分别独立地求解隧道内的气体温度场和围岩温度场 .为克服这一困难，我们利用在洞壁表面上，固体温度等于气体温度这一事实，把隧道内气体的温度和围岩内固体的温度放在一起求解，这样壁面温度将作为末知量被解出来 .只是需要注意两点:解流体温度场时不考虑相变和解固体温度时没有对流项;在洞壁表面上方程系数的光滑化 .另外，带相变的温度场的算法与文献 [3] 相同 .2.3 热参数及初边值的确定热参数的确定方法 : 用 p=1013.25-0.1088H 计算出海拔高度为H 的隧道现场的大气压强，再由P计算出现场空气密度，其中T 为现场大气的年平均绝对温GT度，G 为空气的气体常数 .记定压比热为C P ,导热系数为，空气的动力粘性系数为.按ａ和计算空气的导温系数和运动粘性系数.围岩的热物理C P参数则由现场采样测定 .初边值的确定方法 :洞曰风速取为现场观测的各月平均风速 .取卞导风进曰的相对有效气压为 0，主导风出口的气压则取为p (1 kL / d) v2/ 2[ 5]，这里 k 为隧道内的沿程阻力系数， L 为隧道长度， d 为隧道端面的当量直径，为进口端面轴向平均速度 .进出口气温年变化规律由现场观测资料，用正弦曲线拟合，围岩内计算区域的边界按现场多年冻土下限和地热梯度确定出适当的温度值或温度梯度.3计算实例按以上所述的模型及计算方法，我们对大兴安岭西罗奇 2 号隧道内气温随洞曰外气温变化的规律进行了模拟计算验证，所得结果与实测值 [6] 相比较 ,基本规律一致 .西罗奇 2 号隧道是位十东北嫩林线的一座非多年冻土单线铁路隧道，全长1160 m ,隧道近西北一东南向，高洞口位于西北向，冬季隧道主导风向为西北风 .洞口海拔高度约为 700 m ,月平均最高风速约为 3m/s,最低风速约为 1.7m/s.根据现场观测资料，我们将进出口气温拟合为年平均分别为 -5 0C和-6.4 0C，年变化振幅分别为18.9 0C和 17.6 0C的正弦曲线 .隧道的当量直径为 5.8 m,沿程阻力系数取为 0.025.由于围岩的热物理参数对计算洞内气温的影响远比洞口的风速、压力及气温的影响小得多，我们这里参考使用了大坂山隧道的资料 .图 1 给出了洞口及洞内年平均气温的计算值与观测值比较的情况，从进口到出口，两值之差都小于 0.2 0C .图 2 给出了洞内 (距进出口 l00m 以上 )月平均气温的计算值与观测值比较的情况，可以看出温度变化的基本规律完全一致，造成两值之差的主要原因是洞口气温年变化规律之正弦曲线的拟合误差，特别是1979 年隧道现场月平均最高气温不是在 7 月份，而是在 8 月份 .4对大坂山隧道洞内壁温及围岩冻结状况的分析预测4.1 热参数及初边值按大坂山隧道的高度值 3 800 m 和年平均气温 -3 0C，我们算得空气密度0.774kg / m3；由于大气中含有水汽，我们将空气的定压比热取为[7]C p 1.8744kJ / m s 导热系数 2.0 10 2 W / m 0C ，空气的动力粘性系数取为9.21810 6 kg / m s ,经计算，得出空气的导温系数 a 1.3788 10 5 m2 / s 和运动粘性系数 1.19 10 5 m2 / s .考虑到车体迎风面与隧道端面相比较小、车辆在隧道内行驶速度较慢等因素，我们这里忽略了车辆运行时所形成的活塞效应对气体扩散性能的影响.岩体的导热系数皆按完好致密岩石的情况处理，取岩石的干容重d2400kg / m3时，含水量和末冻水含量分别为W=3% 和W=1 % ，u1.9W / m.o c ,f(0.8 4.128w u )C f1 W2.0W / m.o c 岩石的比热取为C V0.8kJ / kg.0 C ，(0.8 4.128w u )d , C ud .1 W另外，据有关资料，大坂山地区月平均最大风速约为3.5 m/s ，月平均最小大风速约为 2.5m/s;我们将洞口风速拟合为 v t)[0.028 (t7) 2 2.5](m s ，这 ( / )里 t 为月份 .洞内 风 速初 值为 ： U ( 0, x , r ) U a (1( r) 2 ), V ( 0 , x , r ) 0. 这 里取RU a 3.0m / s .而将温度的初边值取为这里记 f (x) 为多年冻土下限到隧道拱顶的距离， Ro = 25m 为求解区域的半径 .地热梯度取为3%，洞外天然年平均气温 A=-3 0 C ，年气温变化振幅 B=12 0C .对于边界 R = Ro ，我们先按第一类边值 (到多年冻土下限的距离乘以 3 %)计算，发现一年后，在半径为 5m 到 25m 范围内围岩的热流方向己经发生转向 .考虑到此后围岩会继续冷却，但在边界R= R 0 上又受地热梯度的作用，我们近似地将边界 R= Ro 作为第二类边界处理，即把由定边值计算一年后 R=R 。