本研究利用COMSOL的传热模块、电流模块、PDE模块以及电流和传热耦合的电磁热模块，构建土壤电阻加热模型。为简化模型概念，在模拟时做如下假设：1)忽略土壤和水在加热过程中密度、热容随温度的变化；2)将土体视为均质且各向同性的多孔介质，土壤初始温度均匀一致；3)忽略土体在加热过程中的热变形；4)忽略电极和土壤之间的接触电阻；5)由于实验土柱较短，需要考虑水分在低于沸点时的挥发^[13]，但在模拟大尺度的场地加热中，忽略水分的挥发；6)忽略水蒸气对热量传递的影响。

电阻加热土壤过程中能量的控制方程如式(1)所示。

式中：ρ_eff为单元的有效密度，kg·m⁻³；C_eff为单元的有效热容，J·(kg·K)⁻¹；T为温度，K；t为时间，s；λ为导热系数，W·(m·K)⁻¹；σ为土壤电导率，S·m⁻¹；E为电势梯度，V·m⁻¹；m_LG为水的气化速率，kg·(m³·s)⁻¹；ΔH_vap为水的潜热，J·kg⁻¹；ρ_f为流体密度，kg·m⁻³；C_f为流体热容，J·(kg·K)⁻¹；u_f为流体流速，m·s⁻¹。

土壤单元的有效密度ρ_eff和有效热容C_eff由土壤中固，液，气三相的体积分数决定，如式(2)~式(3)所示。

式中：θ表示各相的体积分数；ρ为各相的密度，kg·m⁻³；C为各相的热容，J·(kg·K)⁻¹；下标S，L，G表示固，液，气三相。

土壤含水量的控制方程如式(5)所示。

式中：D_L为导水系数，m²·s⁻¹；α为比例常数，m²·s⁻¹；θ_L^*为残余饱和度。

液态水变为气态水的情况可分为2种，一种是低于水的沸点时的挥发，一种是到达沸点时的沸腾，用式(6)可以得到水的气化速率。需要注意，只有水的饱和蒸气压(p*)大于等于外部气压(p_G)且含水量大于0时，沸腾才会发生。可以通过安托因方程^[14](式(7))计算不同温度下水的饱和蒸气压，进而判断温度是否到达沸点。

式中：A、B、C为经验常数；m_vap为挥发速率，kg·(m³·s)⁻¹；k_vap为蒸发速率常数，s⁻¹。

土壤单元的热导率会随着温度和含水量的变化发生极大的变化，在此使用TARNAWSKI等^[15]推导出的经验公式，如式(8)~式(9)所示。

式中:下标sat和dry分别表示饱和和干燥状态的土壤；a~g为经验常数；S_w为水饱和度。

电阻加热一般使用低频率电压(50~60 Hz)，产生的电磁波长远大于系统的物理尺寸，位移电流可以忽略，因此可以假设电阻加热产生的电场为准静态电场，可以通过将欧姆定律代入电流连续性方程中来求解电势分布，电流连续性方程如公式(10)所示。

土壤电导率则使用Archie定律进行计算^[16-17]，考虑温度对电导率的影响^[18]，如式(11)所示。

式中: ψ为电势，V；φ为孔隙度；m、n和β为经验常数(m为胶结系数，n为饱和度系数，β为温度系数)；σ_L为土壤溶液电导率，S·m⁻¹。

模拟单相交流电加热时，其电势分布与直流电基本一致，可以设置2个电极分别为接地和施加的电压。但模拟三相交流电时，电极电势与直流电则完全不同，三相交流电的电极电势可以视为由虚部和实部组成^[19]，如式(12)所示。

式中：ω为角频率，rad·s⁻¹；t为时间，s；θ为相角，°(三相交流电分别为0 °、120 °和240 °)；j为虚部；E₀为正弦交流电振幅的绝对值，一般为电压的2^1/2倍，V。

电阻加热土柱装置模型如图1(a)所示，装置高16 cm，内径4.5 cm。装置的罐体用不导电的聚四氟乙烯制成，内部装填细砂。装置2边为对称分布的用螺丝与螺母固定的不锈钢片电极，螺丝作为导电连接，电线连接到螺丝上以引入电压，热电偶从装置底部插入，用于监测土柱中心点的温度。使用Solidworks对电阻加热土柱装置进行建模，得到电阻加热土柱装置的几何模型，再导入到COMSOL中用于模拟，几何模型的网格划分均采用三角单元网格，单元大小选择细化。

为了模拟实际场地中电阻加热的过程，构建了图1(b)所示的三相电阻加热场地模型。场地模型为长宽15 m、高5 m的长方体，在场地中部按等边三角形放置3根半径20 cm、高5 m的电极，电极间距为6 m，并在3个电极构成的等边三角形中心点设置了温度监测点。

电阻加热土柱装置数值模拟的边界条件为：1)流动边界条件，所有边界均设置为0通量边界条件，这是因为已在方程中添加了描述水蒸发的汇项；2)温度边界条件，所有边界均设置为自然对流热通量边界；3)电势边界条件，2个电极分别设置为接地和电势，其他边界设置为电绝缘边界。土壤的初始温度、含水量、水的电导率和电极电压根据实测得到的初始值设置。对于各项参数的取值如表1所示。

电阻加热场地模拟的边界条件为：1)流动边界条件，模拟地下水位上涨时，底部边界设置为通量边界条件，其他边界设置为0通量边界条件；2)温度边界条件，地下水流入的边界设置为流入边界，流入温度10 ℃，其他边界设置为热绝缘边界；3)电势边界条件，3个电极分别设置为三相电势中的1相，其他边界设置为电绝缘边界。各项参数的设置与土柱装置实验相同。

为验证模型在不同的土壤含水量和土壤溶液电导率时的准确性，使用电阻加热土柱装置进行了验证实验。由于装置较小，电场强度较大，故使用电导率较小的细砂充当模拟土壤。分别向300 g细砂中加入40、50、60、70 mL的去离子水，质量分数为0.3%的NaCl溶液和质量分数为0.5%的NaCl溶液，充分搅拌，并装填进装置中。为了保持砂的孔隙度一致，在装填时保证装填高度为15 cm。加入去离子水的砂在100 V的条件下进行电阻加热，加入NaCl溶液的砂在50 V的条件下进行电阻加热。由于模拟只探究电阻加热的升温过程，所以实验只进行到温度到达100 ℃就停止，未加热到100 ℃时则在加热2 h后停止。

为模拟实际工程环境下土壤内热量传递规律，用实际场地模型探究了地下水流动、电场强度和电极间距的变化对土壤温度变化的影响。模拟的工况如表2所示，其他参数与表1一致，土壤溶液电导率设为0.1 dS·m⁻¹。

图2给出了不同含水率和电导率的情况下，土柱中心热电偶监测的温度变化。可以看出，随着含水量和土壤溶液电导率的增加，中心点加热到100 ℃所需要的时间不断减小。而当含水量较小或者土壤溶液电导率较小时，升温速率较小，甚至出现温度平台的情况。升温速率较小是由于此时土壤电导率较小^[20]，电流产生的焦耳热较小。升温速率减小则是因为，土柱中挥发掉的水份随着加热时间的增长逐渐变多^[21]，进一步降低了土壤电导率，使电流产生的焦耳热进一步减少。此外，由图2(a)可以看到，加入不同体积的去离子水后的细砂依然可以被加热，而去离子水的电导率几乎可以忽略不计。这说明，细砂中的离子溶解进入了去离子水中，提高了去离子水的电导率。

为验证模型的可靠性，将不同含水量和土壤溶液电导率的电阻加热实测值和模拟值进行了对比。利用均方误差MSE和平均相对误差MRE评价模拟值和实测值的差异(式(12)~式(13))。

式中：n为实测数据个数；M_i、S_i分别为第i个实测和模拟得到的数据。

图2(a)中给出的模拟值趋势线与实测值差异较小，实测值和模拟值均方误差为0.05~0.66，平均相对误差为0.42%~1.97%。较小的均方误差和平均相对误差表明，实测数据与模拟数据之间的偏差较小，模型具有较好的准确性。加入不同体积的0.3%NaCl(图2(b))和0.5%NaCl(图2(c))溶液，土壤溶液电导率随着NaCl溶液的体积和质量分数的增加而增大。对比实测和模拟预测结果，均方误差为0.38~12.29，平均相对误差为1.15%~5.32%。这表明，在较宽的土壤电导率范围内(加入70 mL的0.5%NaCl溶液时，土壤的电导率为1.1 dS·m⁻¹)，实测数据与模拟数据之间的偏差也较小，模型用于预测不同电导率的土壤加热过程是可靠的。误差产生的原因可能是，对土壤原有的电导率，实验过程中空气热对流导致的热量散失以及水分挥发速率的估值存在一定的偏差。

电场强度是影响电阻加热升温速率的重要因素，由电极电压和电极间距共同决定，为了考察场地尺度下的温度场以及各条件对温度场的影响，建立了图1(b)的场地模型。图3和图4分别展示了电极间距6 m时，30、60、90 V·m⁻¹电场强度下加热70 d后，温度场的横截面以及位于3个电极中心点(见图1(b))温度的变化。从图3可以看出，电场强度越大，土壤升温速率越快，电极附近的升温更快，温度更高。这与MCGEE等^[10]和HAN等^[22]的研究结果是一致的。电流产生的焦耳热随电场强度的增大而增大，电场强度越大，土壤升温速率越快；电极附近的电流密度最高，产生的焦耳热最多，所以土壤升温速率更快。

从图4可以看到，场地在加热90 d后，30、60和90 V·m⁻¹电场强度获得的中心温度分别为25.5、49.1和100.0 ℃，中心点的平均升温速率分别为0.079，0.420和1.210 ℃·d⁻¹。由式(1)可知，焦耳热与电场强度的平方呈正比，电场强度为30、60、90 V·m⁻¹时中心点升温速率之比应为1∶4∶9，小于模拟得到的比值，即1∶5.3∶15.5。这说明，中心点的温度可能是电流焦耳热和外部热传导叠加共同决定。此外，电场强度为90 V·m⁻¹时，当中心点温度达到97 ℃后，升温速率明显放缓。这是因为，此时的电极温度已经达到水的沸点，电极土壤水分蒸发带走了大量热量，导致土壤热导率和土壤电导率下降，从而使升温速率下降。

电极间距的设置决定了电极井的数量，会极大地影响修复场地的成本和热传导过程。为探究电极间距这一单因素对电阻加热过程的影响，在保持电场强度90 V·m⁻¹不变的条件下，模拟电极间距3、6和9 m时中心温度的变化，结果如图5所示。从图5(a)可以看出,在电场强度为90 V·m⁻¹，电极间距为3、6和9 m时，中心点的温度达到100 ℃的时间分别为84、62和75 d，达到100 ℃的能耗分别为10 418、23 375和51 311 kW·h。综上可知，6 m为最适宜的电极间距，此时中心点的升温速率最快，且相较于电极间距为3 m时场地修复需要布设的电极井数量更少，相较于电极间距为9 m时需要的能耗更少。从图5(a)还可以看出，随着电极间距的增大，中心点升温速率先增大后减小。这是因为，为了保证电场强度一定，电压随电极间距增大而增大，导致电极处的升温速率也随之增大，这一定程度上弥补了间距增大导致的中心点热传导距离增大的不足。但是，9 m间距下的电压增加，造成电极附近温度过早达到水的沸点，导致电极周围土壤水分过早蒸发，土壤热导率和电导率下降，升温速率略为下降。

为了考察电极附近土壤水分蒸发对中心点升温速率的影响，模拟了电场强度90 V·m⁻¹、电极间距为9 m时，对电极附近补水的情况下，中心点温度的变化。从图5(b)可以看到，在加热40 d后，补水的情况下中心点温度明显高于不补水的情况，补水的情况下中心点的温度达到100 ℃的时间为68 d，比不补水的情况下早8 d。此结果表明，电极周围土壤水分过早蒸发是造成升温速率下降的原因之一，补水可以使土壤含水率增大，增加土壤的电导率和热导率，提高了升温速率。这与葛松等^[3]的研究结果一致。监测电极的电流变化，并及时对电极附近进行补水，对于更快加热到目标温度至关重要。

使用原位电阻加热的场地一般修复深度都较大，有必要考虑地下水流动和地下水位变化对温度场的影响。为进行预测，分别设定了0、0.1、0.2、0.3 m·d⁻¹的地下水流速下加热30 d和0、0.05、0.10、0.20 m·d⁻¹的地下水上涨速度下加热30 d的条件进行模拟，流入的地下水出的温度假定为10 ℃，温度场的变化结果如图6和图7所示。从图6可以看到，随着地下水流速的增加，场地左侧的温度逐渐降低到10 ℃，温度场右移的趋势越明显，中心点的最高温度从100.0 ℃降低到38.5 ℃。这说明，较高的地下水流速会将场地中的热量带向下游，对场地升温速率造成不利影响，与MUNHOLLAND等^[23]在二维沙箱中得到的结果一致。实际工程中应尽量减小地下水流速，必要时可增加侧向的物理屏障来形成止水帷幕，或通过蒸汽注入等手段来增加流入地下水的温度。

从图7可以看到，随着地下水位上涨速度提升至0.1 m·d⁻¹，场地下部和两侧的温度逐渐降低到10 ℃，但由不饱和区域变为饱和区域的土壤升温速率上升，中心点最高温度从51 ℃升高到100 ℃。随着地下水位上升速度达到0.2 m·d⁻¹，相较于0.1 m·d⁻¹时，温度场明显向上移动，中心点最高温度从100 ℃降低到54 ℃。这是因为，地下水位上涨可以起到一定的补水作用，提高不饱和区域土壤的电导率，从而提高土壤升温速率，但地下水上涨速度过快还是会导致流失的热量多于通过补水增加的焦耳热，使得场地温度从底部开始下降，不利于场地的修复。应当指出，由于实际场地地层的复杂性，地下水的流速和流向并不是均匀的。因此，对于非均质地层，可以建立相应的多模块的空间耦合模型，通过改变各地层模块的物理参数，从而实现数值模拟目标。

1)电阻加热温度场数值模型与土柱装置实验数据吻合度较好，模拟值和实测值均方误差为0.05~12.29，平均相对误差为0.42%~5.32%，数值模型具有较好的预测准确性。

2)场地模型研究发现，升温速率与电场强度成正比，电极处升温速率最快；电场强度为90 V·m⁻¹时6 m为最适宜的电极间距，此时中心点升温速率最快，能耗相对较低，需要的电极井较少；模型考虑水分蒸发的情况下，电极周围土壤水分会更早蒸发，导致土壤升温速率下降，此时对电极附近补水可以显著增加土壤升温速率，故在实际工程中应监测电流和功率，并及时对电极附近进行补水。

3)场地地下水流动会带走热量，不利于场地的热修复。但对于不饱和场地，地下水位上涨速度小于等于0.1 m·d⁻¹时可以起到补水作用，提高场地升温速率。故对于地下水流速过高的场地应采取水力平衡或设置止水帷幕等工程措施以缓解热量流失。