A 3D discrete FEM iterative algorithm for solving the water pipe cooling problems of massive

See discussions, stats, and author profiles for this publication at: /publication/280874377 A 3D discrete FEM iterative algorithm for solving the water pipe cooling problems of massive concrete structuresARTICLE in INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS · JULY 2015Impact Factor: 1.38 · DOI: 10.1002/nag.2409READS84 AUTHORS, INCLUDING:Zhao LanhaoHohai University32 PUBLICATIONS 15 CITATIONSSEE PROFILEA 3D discrete FEM iterative algorithm for solving the water pipe cooling problems of massive concrete structuresJing Cheng*,†,T.C.Li,Xiaoqing Liu and L.H.Zhao1College of Water Conservancy and Hydropower Engineering,Hohai University,Nanjing,210098,ChinaSUMMARYWater pipe cooling has been widely used for the temperature control and crack prevention of massive con-crete structures such as high dams.Because both under-cooling and over-cooling may reduce the ef ficiency of crack prevention,or even lead to great harm to structures,we need an accurate and robust numerical tool for the prediction of cooling effect.Here,a 3D discrete FEM Iterative Algorithm is introduced,which can simulate the concrete temperature gradient near the pipes,as well as the water temperature rising along the pipes.On the basis of the heat balance between water and concrete,the whole temperature field of the problem can be computed exactly within a few iteration steps.Providing the pipe meshing tool for building the FE model,this algorithm can take account of the water pipe distribution,the variation of water flow,water temperature,and other factors,while the traditional equivalent algorithm based on semi-theoretical solutions can only solve problems with constant water flow and water temperature.The validation and convergence are proved by comparing the simulated results and analytical solutions of two standard second-stage cooling problems.Then,a practical concrete block with different cooling schemes is analyzed and the in fluences of cooling factors are investigated.In the end,detailed guidance for pipe system optimization is provided.Copyright ©2015John Wiley &Sons,Ltd.Received 16May 2014;Revised 15February 2015;Accepted 5June 2015KEY WORDS:numerical methods;iterative algorithm;water pipe cooling;transient heat transfer;thermal stress1.INTRODUCTIONMassive concrete structures such as hydraulic dams,sluices,pumping stations,and nuclear plant foundation always suffer crack problems.Research on crack formation and relevant prevention measures remains an important topic for civil engineers.Despite the fact that dry shrinkage,poor material quality,and many other factors may contribute to the cracking,thermal stress is considered as the toughest to treat with for massive concrete structures.Normally,the cement hydration lasts several days for common concrete or weeks for RCC,during which the core temperature will rise signi ficantly without controlling measures [1,2].However,because of the inherent poor heat conductivity of concrete,it may take several years for core temperature to fall down to steady or quasi-steady state in accordance with natural conditions.During the early age after construction,the core temperature is much higher than the surface temperature,and the inner expansion will cause surface cracks.After years of heat emission,the core temperature drops signi ficantly,and then,the inner shrinkage crack may occur under the constraints of surrounding parts.Apparently,to prevent these two major kinds of thermal cracks,the key issue is to control the temperature rising and the inner-surface temperature difference.Many measures have been taken such as cement reduction,aggregate pre-cooling,placing *Correspondence to:J.Cheng,College of Water Conservancy and Hydropower Engineering,Hohai University,Nanjing,210098,China.†E-mail:mscj042@INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS Int.J.Numer.Anal.Meth.Geomech.(2015)Published online in Wiley Online Library ().DOI:10.1002/nag.2409J.CHENG ET AL.temperature controlling,and block surface protection[3].Among them,the water pipe cooling proves to be the most efficient and economic.Dating from the1930s,the water pipe cooling technique wasfirst used by the Bureau of Reclamation in Hoover Dam,of which the thin-walled metal pipes were embedded in the concrete block during construction period[1].Since then,it spread over the world from common concrete to roller-compacted concrete(RCC), and from dams to sluices and other structures[4–6].In the1970s,a high-density polyethylene(HDPE)pipe with high thermal conductivity arose and soon took over metal pipe in hydraulic projects,such as Sayano–Shushenskaya Dam and Ertan Dam,because of its greatflexibility for shaping[6].Figure1(a)shows the configuration of Jinping-I Arch Dam with26sections divided by temporary transverse joints.Each dam section was placed layer by layer at a vertical interval of3m.The cooling pipes were usually embedded in horizontal section planes like an S curve shown in Figure1(b),with both vertical and horizontal pipe spacing of about1.2–3.0m.To keep the effectiveness,the pipe length should better not exceed250m. For those huge concrete structures like arch dams,which arefirst built by several parts and then joined together to meet their uneven deformations,the cooling process generally includes two stages.Thefirst-stage cooling begins several hours after the concrete placement,aiming to suppress the maximum concrete temperature rising.Before the joint grouting between neighboring parts,a second-stage cooling is carried out to make the inner temperature decrease till it gets near its static state.There are three essential categories of computation methods for pipe cooling problems:the analytical methods,the equivalentfinite element method(FEM),and the discrete FEM.In1949,the US Bureau of Reclamation obtained thefirst analytical solution in series for second-stage cooling problems via the method of separation of variables,with the following assumptions:(1)both the temperature andflow of cooling water were constant;(2)pipes were equally spaced;and(3)exothermal heating was complete and the concrete temperature was evenly distributed at the start of pipe cooling[7].Another solution in series was given out by Zhu in1956using the Laplace transform method,with amuchFEM ITERATIVE ALGORITHM FOR WATER PIPE COOLING PROBLEMS OF CONCRETEhigher convergence rate.After that,Zhu proposed analytical solutions forfirst-stage cooling problems with metal pipes in1957[8]and non-metal pipes in1999[9].Analytical solutions provide quantitative insight into the cooling effect,but it lacks practical utilization.With these analytical solutions,an equivalent FEM was developed for cooling problems,of which the cooling effect was considered as a negative heat source in an average meaning[10,11].This brought great convenience for the application of FEM.Though without obtaining the real temperature gradient and stress near the pipes,the equivalent algorithm remains to be the most widely used method in engineering practice.The discrete FEM directly simulates the pipes with elements;thus,the inner heat sources and surface convection along interior pipe wall could be considered simultaneously.By the early1990s, three methods were used for the computation of water temperature along the pipes,the simple method,the iteration method,and the prediction method[10].However,requirements for small element size near pipes and consequent large element number made it difficult for practical use.In recent years,great progress has been achieved in the controlling of pipe cooling on account of the informationization of construction and the digitalization of dam data.For this new pipe cooling system often denoted as‘intelligent water-pipe cooling’,the real-time information of concrete temperature, water temperature,waterflow,and other indexes are collected automatically,and the cooling process can be adjusted freely[12].Accordingly,the executive cooling criteria extends from a single maximum allowable temperature to a temperature history curve based on material properties and structure features,and the optimization of cooling scheme can bring benefit in millions of dollars for high dams. With the faintness of traditional semi-theoretical FEM for these innovations,many efforts have been made to provide much more accurate and robust numerical tools for the prediction of cooling effect. Xie et ed a3Dfinite element relocating mesh method to simulate the temperature distribution for different pipe cooling schemes of RCC arch dams during construction and operation period[5]. Kim et al.developed a line element to model the pipe,and the internalflow theory was adopted for calculating the temperature variation of pipe water[13].Yazdi et al.developed a temperature simulation model usingfinite volume method[14].On the basis of the heat balance between concrete and cooling water,Liu proposed a direct algorithm for metal pipe cooling problems,of which the fully coupled finite element equations were established with both the piped water temperature and the concrete temperature as unknown.Thus,the water temperature can be solved directly without iteration[15].In2003,Zhu et al.proposed an FEM that uses the discretefinite element to model the pipes,and the water temperature can be calculated by iterations[16].To distinguish this method from other approaches based on FEM,we call it a‘discrete FEM iteration algorithm’.Without any assumptions,almost all factors can be simulated easily,including pipe cooling parameters,concrete material parameters,and complex boundary conditions.However,this method did not attract much attention until recently when the great progress was achieved in the digitalization of dam data and high performance computing.In this paper,the kernel idea of the3D discrete FEM iterative algorithm for pipe cooling problems of massive concrete is presented.The key items,including the discretization of matrix equation,the calculation of water temperature along the pipe,the iterative algorithm for temperaturefield simulation,and the incremental method for thermal stress calculation,are described,respectively.For those practical problems with complex embedded pipes,a very useful and novel pipe meshing tool is developed on the basis of coarse mesh and pipe distribution information.Then,the two standard second-stage pipe cooling problems are simulated for the verification of the algorithm and its convergence evaluation.Finally,a practical application is provided as an instruction for the proposed method and for the practical design of pipe cooling system.2.THEORY AND METHODOLOGY2.1.Discretized matrix equation of heat conductionTo solve the pipe cooling problems of massive concrete structures above,the3D transient heat conduction governing equation is involved.In the given domainΩincluding the structure and its basement within a certain range,it can be expressed as[17,18]∂∂xλx∂T∂xþ∂∂yλy∂T∂yþ∂∂zλz∂T∂zþQÀc cρ∂T∂t¼0(1)where T=T(x,y,z,t)denotes the temperaturefield at time t;λx,λyλz denote the thermal conductivity in the x,y,and z directions,and for homogenous material,λx=λy=λz=λ;The heat source term Q here has the following form for newly placed concrete:Q¼c cρ∂θ∂t(2)whereθ=θ(t)denotes the adiabatic temperature rise of concrete due to hydration heat;This partial differential equation can be solved with the following initial and boundary conditions, which is known or specified:T¼T0x;y;zðÞ;when t¼t0initial condition(3)T¼T b x;y;z;tðÞonΓT;essential boundary(4)Àλ∂T∂n¼βTÀT aðÞonΓq;convection boundary(5)where T b(x,y,z,t)is the known boundary temperature function or specified function,∂T/∂n is the temperature gradient along the n direction,and for3D problems,it can be expressed as follows:∂T ∂n ¼∂T∂x n xþ∂T∂y n yþ∂T∂z n z(6)and n x,n y,and n z are the direction cosines of the external normal to the boundary.Generally,ΓT includes the bottom of rock-base model where the temperature is specified according to geological information,whileΓq includes the surface in contact with water or air.The above problem can be discretized into elements and solved in the following matrix form using standard FEMs[17]:C_dþKd¼P(7)where C and K are the global specific heat matrix and heat conductivity matrix,P is the load vector, and d and_d are the node temperature vector and its time derivative,respectively.C¼∫Ωc cρN T N dΩ(8)K¼∫Ωλ∂N∂xT∂N∂xþ∂N∂yT∂N∂yþ"∂N∂zT∂N∂z#dΩþ∫Γ3βN T N dΓ(9) P¼∫Ωc cρ∂θ∂t NT dΩþ∫Γ2N T q dΓþ∫Γ3N TβT a dΓ(10)d¼T1T2…T nf g T n is the total node number(11) N is the shape function described in the FEM method.J.CHENG ET AL.Using the weighted residual approximation for temperature field within the time interval Δt =t n+1Àt n ,Eq.7can be written asC =Δt þK ϕðÞd n þ1þÀC =Δt þK 1ÀϕðÞ½ d n ¼P (12)This is the well-known θalgorithm or θmethod for solving the first-order equations.To avoid confusion with the previous mentioned heat generation,here,the weighting parameter ϕis used instead.Different values of φlead to various time difference approximations:(1)ϕ=0,forward-difference approximation (or Euler difference approximation);(2)ϕ=1/2,central-difference approximation (Crank –Nicholson difference approximation);(3)ϕ=2/3,Galerkin difference approximation;;and (4)ϕ=1,backward-difference approximation.Then,the solution of Eq.12yields the following:C =Δt þK ϕðÞd n þ1¼P ÀÀC =Δt þK 1ÀϕðÞ½ d n ¼P þC =Δt þK ϕÀ1ðÞ½ d n(13)It can be written as follows:K C 1d n þ1¼P þK C 2d n(14)where K C1and K C2are working matrices:K C 1¼C =Δt þK ϕ;K C 2¼C =Δt þK ϕÀ1ðÞ(15)The load vector can be calculated as follows:P ¼P n 1ÀϕðÞþP n þ1ϕ(16)Here,ϕ=1is retained in the simulation of this paper.2.2.Calculation of water temperature along the pipe`On the early stage of arti ficial pipe cooling,the temperature difference between concrete and water is comparatively large,which leads to the signi ficant rising of water temperature along the pipes.The water temperature in the inlet of the pipe is preset from refrigeration plant or just from the river because of its source;thus,for solving the pipe cooling problems de fined by Eq.1under conditions Eqs.3–5,the only left to be calculated is the water temperature increase along the pipe.For a typical vertical section of the concrete block in Figure 2(a),the pipes are distributed at a spacing of S H ×S V ,and each pipe only effects a certain part approximating to a column or prism in Figure 2(b,c).As a type of convection boundary de fined in Eq.5,the heat transferfrom Figure 2.Physical model for concrete block embedded with cooling pipe.FEM ITERATIVE ALGORITHM FOR WATER PIPE COOLING PROBLEMS OF CONCRETEconcrete to water can be denoted as follows:q¼Àλ∂T∂n(17)The water temperature increase between two neighboring sections S1and S2of the column is denoted asΔT S.During a very small time interval d t,1.The heat transferred from concrete to water through the inner surfaceΓp of the pipe betweenSection S1and S2can be expressed as follows:d Q p¼∬Γp q i d s d t¼Àλ∬Γp∂Td s d t(18)2.The heat energy of the waterflowing through the inlet section S1isd Q S1¼c wρw q w T S1d t(19)3.The heat energy of the waterflowing through the outlet section S2isd Q S2¼c wρw q w T S2d t(20)In Eqs.19and20,T S1and T S2are the water temperatures in Sections S1and S2.4.The change of the heat energy of the water between S1and S2is as follows:d Q w¼∫S2S1c wρw∂T S∂t d tA p d l(21)where A p is the area of the pipe section,T s denotes the water temperature in any section S,and d l is the distance between S1and S2.According to the heat equilibrium,we haved Q S2¼d Q S1þd Q PÀd Q w(22)Substituting Eqs.18–21into Eq.22,the water temperature increase between S1and S2can be expressed as follows:ΔT S¼T S2ÀT S1¼Àλc wρw q w∬Γp∂T∂n d s d tÀA pq w∫S2S1∂T S∂t d td l(23)Considering that both the volume and temperature increase of water are very small within d l, d Q w is a higher order trace relative to d Q p and can be omitted;thus,Eq.23can be written as follows:ΔT S¼Àλc wρw q w∬Γp∂T∂n d s d t(24)Providing the concrete temperaturefield is known,the above integral can be executed easily on Γp.Because the inlet water temperature T w0is known,we can calculate the water temperature inJ.CHENG ET AL.any section S i by dividing the pipe into m sections and accumulating the temperature increase along the pipe as follows:T S i ¼T w0þ∑i j ¼1ΔT S j ;i ¼1;2;3;…;m (25)2.3.Iterative algorithm for temperature fieldFrom Eq.24,it can be seen that the increase of water temperature along the pipe is in fluenced by the axial temperature gradient of the concrete.Note that the water temperature and the structural concrete temperature interact with each other;we adopt the following iterative algorithm to solve this nonlinear boundary problem.1.At the first step,the water temperature along the pipe at time t =t i is assumed to be constant as equal to the inlet water temperature T w0;then,a temporary transient temperature field T 0(x ,y ,z ,t i )can be calculated.2.On the basis of the computed results T 0,the radial concrete temperature gradient along the pipe (∂T /∂n)0is obtained,as well as the water temperature along the pipe T S 0using Eqs.24and 25.3.With the updated water temperature T S 0,the transient temperature T 1(x ,y ,z ,t i )can be calculated,which is much more close to the real temperature than T (x ,y ,z ,t i ).4.By repeating steps (2)and (3),we can obtain a series of (T S k ,T k (x ,y ,z ,t i ),k =0,1,2,…).When two neighboring series are close enough to ful fill the following condition,stop the iteration,max T k þ1Si ÀT k SiÀÁ<ε;i ¼1;2;3;…;m ;ε>0(26)where k is the iteration times and εis the speci fied tolerance,which can be set as 0.01°C.Normally,2–3iterations will be enough.2.4.Pipe meshing toolThe meshing for discrete FEM algorithm is very time consuming.Here,an effective meshing tool PipeMeshF is developed to create the pipe-embedded mesh model on the basis of a coarse FEM model and the pipe distribution information including inlet elements,corner elements,and outlet paring a coarse model and its associated pipe-embedded model in Figure 3,we can observe the following:1.With the PipeMeshF tool,only those elements that pipes go through are remeshed;thus,the total element number is not enlarged somuch.Figure 3.Finite element meshing for pipe-embedded model using PipeMeshF .FEM ITERATIVE ALGORITHM FOR WATER PIPE COOLING PROBLEMS OF CONCRETE2.Thefine mesh model with embedded pipes can be easily created from a coarse model by any pop-ularfinite element meshing packages,provided that those elements pipes go through are hexahe-dron elements.3.Because most parts of the two models keep the same,the results from discrete algorithm withpipe embedded model and equivalent algorithm with coarse model can be compared conveniently.4.The element density around the pipes can be easily changed,depending on the accuracyrequirement.2.5.Thermal stress calculation for massive concreteBecause thefinal purpose of cooling is to control the stress of massive concrete,here,we briefly introduce the theory used for thermal stress calculation.On the basis of the common FEM approach for mechanical problems,the main idea of time-dependent thermal stress calculation is executed by introducing the following relationship between stress and temperature strain:Δσf g¼DΔεÀΔεTf g(27) where{Δσ}is the stress increment,{Δε}the total strain increment,{ΔεT}the total strain increment, and D the elasticity matrix.Concrete creep and time-dependent elastic modulus are the two fundamental features to distinguish thermal stress of concrete structures from that of others.The modulus of elasticity can be expressed by one of the following two equations:EτðÞ¼E01ÀeÀAτBEτðÞ¼E0τHþτ(28)and the unit creep(creep produced under the action of unit stress)Cr(t,τ)is expressed asCr t;τðÞ¼∑n ci¼1ψi1ÀeÀr i tÀτðÞh i(29)ψi¼f iþg iτÀp i(30) where E0,A,B,H,f i,g i,p i,r i,and n c in Eqs.28–30are material constants.With these two factors considered,we call the thermal stress as elasto-creeping or viscoelastic.To reduce the large memory requirement for recording the stress history in the process of creep strain computation,B.F.Zhu proposed an implicit incremental method with unequal time intervals as shown in Figure 4.For this method,the relation between stress increment vector andstrainFigure4.Incremental method for stress calculation.J.CHENG ET AL.increment vector for complex stress state is derived:Δσn f g ¼D n ÂÃΔεn f g Àηn f g ÀΔεT nÈÉÀÁ(31)where D n ÂüE n Q ½ À1;E n ¼E τn ðÞ1þE τn ðÞC t n ;;τn ðÞ;τn ¼τn À1þτn ðÞ=2(32)and ηn is an intermediate variable for calculation the creep strain:ηn ¼∑n c i 1Àe Àr i Δτn ÀÁωi ;n ÈÉωi ;n Èɼωi ;n À1ÈÉe Àr i Δτn À1þQ ½ Δσn À1f g ψi τn À1ðÞe À0:5r i Δτn À1ωi ;1ÈɼQ ½ Δσ0f g ψi τ0ðÞ(33)Here,only the main formulas are collected to save space;more detail about matrix [D],[Q],[Q]À1and the deduction of the implicit method can be found in Reference [6].3.METHOD VERIFICATION AND DISCUSSIONAs indicated in Section 2.2,the concrete region controlled by each pipe can be approximated to a column or prism.Here,a second-stage arti ficial cooling problem for a single concrete column is analyzed by the proposed method to investigate its performance.The problem ful fills the three assumptions mentioned in the third paragraph of introduction,and the column section is shown in Figure 5(a).With an isolated external boundary and without heat source,it possesses analytical solutions for both conditions of metal pipes and HDPE pipes.The initial temperature difference between the concrete and cooling water is de fined as T 0.Other parameters are taken as c =0.016and b /c =100.For the numerical implementation based on the proposed method,four schemes (named as M1T1,M1T2,M2T1,and M2T2,respectively)with different mesh models (named as Mesh1and Mesh2)and time marching schemes (named as TList1and TList2)are computed and compared.Both mesh models are created in a way that the element size d increases along the axis with a constant ratio.As shown in Figure 5(b,c),N Div means the division number of the radius,d max the maximum element size along the axis,and d min the maximum element size along the axis.Two time marching schemes are set asfollows:Figure 5.A single concrete column embedded with cooling pipe.(a)Schematic of section plane;(b)Mesh1,N Div =20,d max /d min =10,d min =0.0198m;(c)Mesh2,N Div =40,d max /d min =20,d min =0.00614m.FEM ITERATIVE ALGORITHM FOR WATER PIPE COOLING PROBLEMS OF CONCRETETList1:0.05,0.15,0.25,0.5,1,1.5,2,2.5,3,4,5,6,8,10,12,14,17,20,25,30,35,40,50,60,70,80,90,100;(unit:day)TList2:0.01,0.02,0.03,0.04,0.05,0.075,0.1,0.125,0.15,0.175,0.2,0.225,0.25,0.3,0.35,0.4,0.45,0.5,0.75,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3,3.5,4,4.5,5,5.5,6,7,8,9,10,11,12,13,14,16,18,20,24,28,32,36,40,45,50,55,60,65,70,75,80,85,90,95,100;(unit:day)Thus,M1T1means the combined scheme with Mesh1and TList1,and similarly M1T2with Mesh1and TList2.The numerical solutions are compared and analyzed for metal pipes and HDPE pipes,respectively.3.1.Second-stage arti ficial cooling with metal pipesWhen the pipe in Figure 5(a)is made of metal,the internal boundary for the concrete column is almost equal to the water temperature because of the high conductivity coef ficient of metal.Analytical solution of temperature field is given by [6]via the Laplace transform method:T r ;τðÞ¼T 0∑∞n ¼12e Àα2n b 2a τ=b 2αn b ÁJ 1αn b ðÞY 0αn r ðÞÀY 1αn b ðÞJ 0αn r ðÞR αn b ðÞ(34)Figure 7.Analytical temperature solution along the radius at different times (N order =5,e T r =2.64e À2).Figure 6.Analytical temperature solution along the radius at different times (N order =3;e T r =3.91e À2).R αn b ðÞ¼c =b ðÞJ 1αn b ðÞY 1αn c ðÞÀJ 1αn c ðÞY 1αn b ðÞ½þJ 0αn c ðÞY 0αn b ðÞÀJ 0αn b ðÞY 0αn c ðÞ½ (35)where both the temperature field T (r ,τ)at t =τand initial temperature T 0at t =0are de fined on a relative temperature coordinate with water temperature as the origin;J 0and J 1are the first and second orders of the Bessel function of the first kind;Y 0and Y 1are the first and second orders of the Bessel function of the second kind;and αn b is the root of the following characteristic equation:J 1αn b ðÞY 0αn c ðÞÀJ 0αn c ðÞY 1αn b ðÞ¼0(36)Note that the analytical solutions in Eq.34are given in series and that the different numbers of order N order we retain may greatly in fluence their accuracy.The solutions with N order =3,5,and 10for this problem are shown in Figures 6–9where a normalized temperature T /T 0,distance r x /b and time a τ/D 2are used for ing the following relative error notation,e Tr ¼T =T 0ðÞmax À1(37)we obtain the relative error of analytical solutions with N order =1,2,3,4,5,10,and 20as e Tr =5.51e À2,4.19e À2,3.91e À2,3.29e À2,2.64e À2,1.66e À3,and 9.3e À7.It can be seen that N order has great effecton Figure 9.Analytical temperature history for different points along the radius.(N order =10).Figure 8.Analytical temperature solution along the radius at different times (N order =10;e T r =1.66e À3).the accuracy of the temperature especially when the time a τ/D 2is small.Because N order >10,the relative error e Tr is less than 0.167%,which can be enough for engineering purposes.Here analytical solutions with N order =20are taken for comparison with the numerical results in Sections3.1and 3.2.The temperature results at three time points for different computation schemes are compared in Figures 10–12,as well as the relative error,which is de fined as follows:e T ¼T =T 0ðÞsim ÀT =T 0ðÞana (38)where the subscript ‘sim ’means simulation results and ‘ana ’means analytical solutions with N order =20.We can conclude the following:(1)simulated results for all schemes are close to analytical solutions,which proves the ef ficiency of the proposed algorithm;(2)for a certain time step list,when the mesh is re fined,the results become much more accurate and closer to the analytical solutions;and (3)for a certain mesh model,when the time step is decreased,the results become much farther from the analytical solutions.This phenomenon of convergence seems strange but can be clari fied well as follows.Remark 1Errors of numerical methods generally include two parts:one due to mesh model and the other from time marching schemes.Here,for the second-stage cooling problem and for the backward-difference approximation (ϕ=1),Eq.12can be deduced asFigure 10.Temperature along the radius and its relative error.(a τ/D 2=0.0781).Figure 11.Temperature along the radius and its relative error (a τ/D 2=0.3906).。