第10期 周扬等：考虑未冻水单管冻结温度场解析解 1651 2问题的求解 Th-V -Ap-Ao）+ 在土体的每一相中，如下形式的温度分布满足基 △0 本方程式(2)： k exp()(15) (i=0,1,2,…,0 式(13)~(15)构成参数入的N阶封闭非线性 方程组，利用牛顿迭代法求解这些方程组后得到A, 6 将,代入式(8)，(10 (12)中获得S,(),a,b,再 其中，4：b,为待定系数，并且有指数积分函数： 将这些参数代入式(6)便可以得到各相中的温度解 Ei()=f 答。 利用初始条件及边界条件式(3)，可以得到 3算例分析 3.1参数的计算 土体孔隙率为0.375，初始体积含水量 由于在两相交界面上任意时刻均需要满足式 为0.2208，土体的冻结特性按式(1)计算，其中S= (4),(5),因此移动边界必须具有如下形式： 60m1.0.=1670lg/m3 S(t)=2λ(a-t)1 (i=1,2,…,N0(8) 各相士体的导热系数按照其成分进行计算，采用 其中，入，为待定系数。 如下公式 利用相界面上的温度关系式(4)，可以得到 k=经经 (16) -1+6-Ei(-A）=T 式中，k为冰的导热系数，2.32W/(m·K):a为冰的 a +bEi(-)= (9) 体积含量：k。为水的导热系数，0.58W1(m·K):k.为 (i=1,2,…,N) 土骨架的导热系数，0.907W1(m·K):日为土骨架体 式中，u,=(ae1/a)n。 积含量：k,为空气的导热系数，0.029W1(m·K):0 利用式(9)，可以将a,b,用A,表示，即 为空气的体积含量 4,=4kT-0E(-A 各相土体体积热容采用式(17)四计算： C.=CB.+C+C+C (17) 4Tko (10) 其中，C为体积热容：C,C,C,C分别为冰，未冻 0 b=4 水，士骨架，空气的体积热容，分别取为1.928x10 a,-TEi(-A)-TEi(-A 4.18×10°，2.2279×10,15.5818×102W1(m3·K)。 若无说明，下文算例中所用参数均与本节相同。 2(λ） (11) 3.2解析解的验证 (i=1,2,…,N-1) 为了验证将十体连相变程离化后得到的 温度场解析解的正确性，将该解析解 土体连续相 模型下的数值解答进行了对比，柱坐标系中土体连续 T-V (12) ,i(-A】 相变过程的控制方程为 式中，△a=T-T:2(a,)=Ei(-A)-Ei(-Aiw)。 c,=-p. (18) 利用相界面上的潜热关系式(5)及的形式可 [0.=f(T) (T≤T) 以 (19) 18.=8 (T>T) -与说p(A》m,+是× △8 式中，k为导热系数 初始、边界条件及计算参数与相应解析解相同 exp(-A)=L'Aiad (13) 对于这样的定解问题，没有解析解，作为对比参考，采 用Matlab中pdepe命令求解其数值解。 下面给出两组算例中解析解与数值解的对比，算 exp(-A)=LA&-1(i=2,3,…,N-1) 例1中0=50W/m,V=2℃，S=60m1:算例2中0 100W/m,V=4℃，=30m. 994-016 China Academic Joumal Electronic Publishing House.All rights reserved.http://www.cnkine第 10 期 周 扬等: 考虑未冻水单管冻结温度场解析解 2 问题的求解 在土体的每一相中，如下形式的温度分布满足基 本方程式( 2) : Ti ( r，t) = ai + biEi － r 2 4αi ( )t ( i = 0，1，2，…，N) ( 6) 其中，ai，bi 为待定系数，并且有指数积分函数: Ei( z) = ∫ z －" eu u du 利用初始条件及边界条件式( 3) ，可以得到 aN = V，b0 = Q 4πk0 ( 7) 由于在两相交界面上任意时刻均需要满足式 ( 4) ，( 5) ，因此移动边界必须具有如下形式: Si ( t) = 2λi ( αi－1 t) 1 /2 ( i = 1，2，…，N) ( 8) 其中，λi 为待定系数。 利用相界面上的温度关系式( 4) ，可以得到 ai－1 + bi－1Ei( － λ2 i ) = T* i ai + biEi( － λ2 i ω2 i ) = T* i ( i = 1，2，…，N) { ( 9) 式中，ωi = ( αi－1 /αi ) 1 /2 。 利用式( 9) ，可以将 ai，bi 用 λi 表示，即 a0 = 4πk0T* 1 － QEi( － λ2 1 ) 4πk0 b0 = Q 4πk0      ( 10) ai = T* i Ei( － λ2 i+1 ) － T* i+1Ei( － λ2 i ω2 i ) Ω( λi ) bi = Δθi Ω( λi ) ( i = 1，2，…，N － 1)        ( 11) aN = V bN = T* N － V Ei( － λ2 Nω2 N) { ( 12) 式中，Δθi = T* i+1－T* i ; Ω( λi ) = Ei( －λ2 i+1 ) －Ei( －λ2 i ω2 i ) 。 利用相界面上的潜热关系式( 5) 及 bi 的形式可 以得到 － k1 Δθ1 Ω( λ1 ) exp( － λ2 1ω2 1 ) ω1 + Q 4π × exp( － λ2 1 ) = L* λ2 1α0 ( 13) － ki Δθi Ω( λi ) exp( － λ2 i ω2 i ) + ki－1 Δθi－1 Ω( λi－1 ) × exp( － λ2 i ) = L* λ2 i αi－1 ( i = 2，3，…，N － 1) ( 14) － kN T* N － V Ei( － λ2 Nω2 N) exp( － λ2 Nω2 N) + kN－1 ΔθN－1 Ω( λN－1 ) exp( － λ2 N) = L* λ2 NαN－1 ( 15) 式( 13) ～ ( 15) 构成参数 λi 的 N 阶封闭非线性 方程组，利用牛顿迭代法求解这些方程组后得到 λi， 将 λi 代入式( 8) ，( 10) ～ ( 12) 中获得 Si ( t) ，ai，bi，再 将这些参数代入式( 6) 便可以得到各相中的温度解 答。 3 算例分析 3. 1 参数的计算 土 体 孔 隙 率 为 0. 375，初始体积含水量 为 0. 220 8，土体的冻结特性按式( 1) 计算，其中 S = 60 m－1 ，ρd = 1 670 kg /m3 。 各相土体的导热系数按照其成分进行计算，采用 如下公式［15］: k = kθi i kθu u kθs s kθa a ( 16) 式中，ki 为冰的导热系数，2. 32 W/( m·K) ; θi 为冰的 体积含量; ku 为水的导热系数，0. 58 W/( m·K) ; ks 为 土骨架的导热系数，0. 90 7 W/( m·K) ; θs 为土骨架体 积含量; ka 为空气的导热系数，0. 029 W/( m·K) ; θa 为空气的体积含量。 各相土体体积热容采用式( 17) ［15］计算: CV = Ciθi + Cu θu + Csθs + Caθa ( 17) 其中，CV 为体积热容; Ci，Cu，Cs，Ca 分别为冰，未冻 水，土骨架，空气的体积热容，分别取为 1. 922 8×106 ， 4. 18×106 ，2. 227 9×106 ，15. 581 8×102 W/( m3 ·K) 。 若无说明，下文算例中所用参数均与本节相同。 3. 2 解析解的验证 为了验证将土体连续相变过程离散化后得到的 温度场解析解的正确性，将该解析解与土体连续相变 模型下的数值解答进行了对比，柱坐标系中土体连续 相变过程的控制方程为 CV T t = 1 r  r kr T ( ) r － Lfρw θu t ( 18) θu = f( T) ( T ≤ Tf ) θu = θ0 ( T ＞ Tf { ) ( 19) 式中，k 为导热系数。 初始、边界条件及计算参数与相应解析解相同， 对于这样的定解问题，没有解析解，作为对比参考，采 用 Matlab 中 pdepe 命令求解其数值解。 下面给出两组算例中解析解与数值解的对比，算 例 1 中 Q= 50 W/m，V = 2 ℃，S = 60 m－1 ; 算例 2 中 Q= 100 W/m，V = 4 ℃，S = 30 m－1 。 1651