基于有限差分法的冻结过程数值分析

基于有限差分法的冻结过程数值分析

在中国东北部和西北部的寒冷地区，土木工程方兴不断，人工土壤冷冻技术在土木工程施工中发挥着越来越重要的作用。要分析自然冻土和人工冻土的冰冻工程，必须正确规定土壤温度变化规律。例如，konrad分凝冰模型和daiyuwa的冻胀模型在计算冻胀前必须正确计算冻面的温度场。然而，在使用有限差分法求解温度场的过程中，由于时间间隔不正确，很容易导致“异常日期”。在这项工作中，我们提出了一种基于自适应时间间隔、半隐藏和全隐藏格式的交替使用方法。这种方法可以有效。防止传统计算过程中出现的相变遗漏,并且计算速度较之固定时间步长方法有了很大的提高.1显热容法求解常见的一维土壤冻结过程热传导方程为,(1)式中:CV,λ为土壤体积热容量(J/(m3·K))和热导率(J/(m·s·K));L,T为融化潜热(J/m3)和温度(℃);ρi,θi为冰密度(kg/m3)和冰含量(m3/m3);z,t为空间坐标(m)和时间(s).求解上述方程的一种有效方法是显热容法,即在发生相变的一个小的温度范围之内构造比热函数,以代替以上方程的右端第二项,只以温度为待求函数,(2)其中:式中:C*V为等效热容;λ*为等效导热系数;Cf,λf分别为已冻土的体积热容和导热系数;Cu,λu分别为未冻土的体积热容和导热系数;Tf-ΔT≤T≤Tf+ΔT为假设小的相变温度范围.2相变影响下的半隐格式方程(2)是与时间和空间都有关系的发展方程,采用有限差分法能对其有效求解.在常用的有限差分格式中,显格式精度最差,隐格式次之,而Crank-Nicolson格式精度较高;但是从解的稳定性角度分析,则隐格式的稳定性最好,因此本文采用隐格式进行求解.根据方程系数的取法不同,隐格式又有全隐格式和半隐格式之分.全隐格式公式如下式中j,n分别为空间和时间坐标标记.所谓半隐格式是将式(3)中系数C*v和λ*的时间步n+1换成n,通过半隐格式只需求解线性方程组即可.事实上,半隐格式和全隐格式对本方程的求解精度影响不大可忽略不计.作者也分别编制了两种不同格式的求解程序进行了验证.由式(2),(3)可见,相变的影响在模拟中是通过在假设相变范围的等效热容来实现的,因此要全面考虑相变的影响,必须保证在每一计算时间步内所有发生相变的空间单元的等效热容能够包含相变潜热项.而土体空间单元的等效热容是否应该包含相变潜热项,则是通过该单元的温度值来判断的.由式(2)知,只有其温度落在Tf-ΔT≤T≤Tf+ΔT范围之内,该单元的相变影响才能被考虑.这样,就不免会出现‘相变遗漏’的情况.比如空间某单元在某一时间步之内温度由T1>Tf+ΔT变为T2<Tf-ΔT,显然该单元在这一时间步内发生了相变,但是根据T1或T2都无法判断出该单元发生相变这一事实,因此单元等效热容也就没能包含本应包含大相变潜热项.由于半隐格式是根据时间步初的温度值、全隐格式是根据时间步末的温度值来判断等效热容是否包含相变潜热项的,因此对于半隐和全隐格式,一旦相应的温度判据即时间步初的温度值和时间步末的温度值没能落在Tf-ΔT≤T≤Tf+ΔT范围之内,该单元的相变影响才能被考虑,否则就会出现下述的‘相变遗漏’现象.例如在任意一个计算时间步内,假设计算空间上某一点在时间步初的温度值为T1,在时间步末的温度值为T2.如果时间步长控制不好,不同的差分格式就会有不同的‘相变遗漏’.采用半隐格式,若T1>Tf+ΔT而Tf-ΔT≤T2≤Tf+ΔT(时段初未冻而时段末处于过渡期)则相变作用就会遗漏;采用全隐格式,若Tf-ΔT<T1<Tf+ΔT而T2<Tf-ΔT相变作用也会遗漏.但是如果T1>Tf+ΔT且T2<Tf-ΔT则不论采用半隐或是全隐格式相变作用都会被遗漏掉.可见如果仅仅采用一种差分格式,那么无论是半隐格式还是全隐格式,对于本文所述问题都会不可避免地存在‘相变遗漏’这一事实,因此是不正确和不可取的.由于冻结过程伴随着大量相变潜热的释放,所以原来采用固定格式算法得出温度将会偏低,而冻结锋面推进速度则会偏快.3速度速的数值分析步骤为了避免数值分析中相变作用的遗漏,必须有效地控制时间步长.本文根据半隐格式和全隐格式的不同对方程的求解精度影响不大这一特点,提出了自调节时间步长、半隐和全隐格式交替使用的方法(ModifiedFDM以下简称MFDM),其中半隐格式直接求解线性方程组即可,全隐格式则要采用迭代法求解.由于本文提出的方法采用的是稳定的计算格式,而且采用了自调节时间步长的技术,所以该方法具有稳定和计算速度快的特点.采用该方法进行数值分析的步骤如下:1)对第一时间步的计算进行特殊处理:第一时间步采用较小的步长,由于第一时间步内温度变化比较剧烈,而且土体从初始的正温开始变化,所以先采用半隐格式进行计算,然后令所有温度落在T<Tf+ΔT的空间点的温度为Tf,再对第一时间步进行全隐格式计算;2)将上一时间步末的空间点按其温度值分为两类:一类是温度落在T≤Tf+ΔT的空间点,另一类是温度大于Tf+ΔT的点;同时采用与上一时间步相同步长的时间步,用半隐格式对本时间步进行试算;3)对第2)步中划分出的温度大于Tf+ΔT的点的试算结果进行判断:如果这些点中的温度最小值(在自上而下冻结的情况下,此点应为这些点中最上端的点)小于Tf-ΔT,就减小时间步长,再用半隐格式对此时间步进行试算;若此最小值大于Tf+ΔT,则增大时间步长,再用半隐格式对该时间步进行试算,直至该最小值落在[Tf-ΔT,Tf+ΔT]之中;4)用第3)步中确定下的最终时间步长对2)中划分出的温度大于Tf+ΔT的点进行全隐格式计算;5)重复2)～4)步,直至计算时间结束.可以设想,采用了上述改进办法,相变遗漏避免了,对于冻结问题也就是说释放出来的相变潜热被全面地考虑了.采用固定格式算法得出的温度值和冻结锋面都会得到相应的纠正,温度会相应升高,锋面推进会相应减缓.4时间步初的温度值采用文献提供的基于FEM的经典算例.空间总长为1m,计算总时间为2.88×105s,采用第一类边界条件,TC=-20℃,TW=10℃;初始温度T0=10℃.参数如下:Tf=0℃,ΔT=1℃,L=338000000(J/m3),λf=2.22(J/(m·s·℃)),λu=0.556(J/(m·s·℃)),Cf=1762000(J/(m3·℃)),Cu=4226000(J/(m3·℃)).1)采用固定时间步长半隐差分格式(FDM)求解.将空间等步划分为400个单元,时间等步划分为192个时间步长.计算结果与文献提供的结果几乎相吻合.通过比较可以看出多处存在相变遗漏,其中计算时间初始阶段出现较多,随着计算时间的增加相变遗漏现象逐渐变少,见表1和表2.表1中在1.5×103～3×103s时间步之内,0.02～0.0275m处4个空间点显然发生了相变,且它们在时间步初的温度值都大于,也即超出了考虑相变的温度范围,由于采用半隐格式,根据时间步初的温度值判断的结论是这4个空间点都没有发生相变(表中黑体数字),因此它们的等效热容将不会包含相变潜热项,也就是说,本时间步之内这4个点的相变都被遗漏掉了.同样在3×103～4.5×103s时间步之内,0.03m和0.0325m两点的相变被遗漏掉了;在4.5×103～6×103s时间步之内,0.035m和0.0375m两点的相变被遗漏掉了;在6×103～7.5×103s时间步之内,点0.04m处的相变被遗漏掉了.在表2中因为相变不像表1所示的时间初那样剧烈,因此相变遗漏相应减少,只有黑体数字表示的两点发生相变遗漏.2)用MFDM求解将空间变步划分为113个单元,初始时间步长设为0.25s,采用此种方法只计算了68个时间步计算就结束了,计算结果见表3.由于对每一时间步内黑体数据以上的点采用半隐格式,而黑体数据及其以下的点则采用全隐格式,所以‘相变遗漏’被有效避免了.比如在表3中的3.17×103～3.68×103s时间步之内,0.025～0.028m处3个空间点采用了半隐格式,点0.03m处则采用了全隐格式,它们相应的温度判据都落在了[Tf-ΔT,Tf+ΔT],因此实际发生的相变都被考虑到了.在3.17×103～3.68×103s时间步之内,点0.025m处已不存在相变,点0.026m和0.028m处均存在相变,对这两个点必须用半隐格式才能考虑相变的影响,而点0.03m和0.031m处因为时间步初末的温度值都在[Tf-ΔT,Tf+ΔT],采用半隐或全隐格式都能考虑到相变影响,点0.033m处也存在相变但只有时间步末的温度值在相变范围内,因此只能采用全隐格式.根据设计的算法在0.033m点以上各点采用半隐格式,在0.033m点(包括该点)以下各点采用全隐格式,是可以满足全部考虑相变影响要求的.图1(横轴表示距土体低温端的空间距离,纵轴表示温度值)比较了用FEM,FDM和本文方法计算所得结果中的两个时刻的温度场.可以看出,FEM和FDM计算结果几乎相同,而用本文方法计算的土体中同一点的温度值则比FEM和FDM的均要高.这符合本文前述的分析和设想:传统FEM和传统的固定格式的FDM均存在相变遗漏,致使没能全面考虑释放的相变潜热,因此计算出的温度偏低.从图1还可看出传统算法和修正算法之间最大的误差出现在冻结锋面附近,这也是符合实际情况的.因为实际冻结过程中相变发生最剧烈的地方就是冻结锋面附近,固定格式算法遗漏相变最多的地方自然也在此处.由于本文提出的方法考虑了原来固定时间步长法没有考虑的‘相变遗漏’的影响,所以比传统方法更加全面地考虑了水相变成冰释放出来的潜热,这必然延缓冻结锋面的推进速度.表4比较了采用半隐格式固定时间步长方法和采用自调节时间步长方法求出的冻结锋面的变化情况.从图2(横轴表示冻结时间,