一种地下硐库热力学与硐壁热传导耦合计算方法

技术领域

本发明专利属于储能领域，涉及压缩空气储能电站地下人工储气硐库充放气过程中的热力学过程计算，以及硐库内空气温度变化引起的硐壁热传导过程，具体为一种地下硐库热力学与硐壁热传导耦合计算方法。

背景技术

压缩空气储能技术是一种大型、长时间储能发电技术，其技术原理如图3所示。在用弃风、弃光及低谷电阶段，压缩空气储能电站利用电能驱动压缩机压缩空气，将其送往地下储气库中储存；在用电高峰期到来时，储气库中的高压空气经过换热器或燃烧室加热，被送至膨胀机内膨胀做功，驱动发电机发电，从而实现对电力能源的削峰填谷等功能。

商业化的压缩空气储能电站地下储气库目前主要为盐穴储气库。它是利用水溶开采方式在地下较厚的盐层或盐丘中注入淡水，将盐岩层溶解，再将饱和或近饱和卤水排出，形成的地下洞穴。盐穴储气库具有密封性好、稳定性高的优点。但盐穴储气库受地理环境限制明显，需建设压缩空气储能电站的地区可能无盐岩分布，无法建造盐穴储气库。

除了盐穴储气库，近年来地下人工硐库也是压缩空气储能电站的一个重要选择。它是在地下硬质岩中人工开挖形成的具有一定体积的储气构筑物。根据气库形状，可以分为大罐式和隧道式两种(见图4～图5)。

不管是何种形式的地下人工硐库，其硐壁结构从硐库内向外主要由密封层、混凝土衬砌和围岩三部分组成。这种形式的储气库气密性由密封材料，如钢材提供，硐库稳定性和变形控制由混凝土衬砌提供。

人工硐库的优势主要体现在工程可控性强，可根据建设需求灵活建造。地下人工硐库是现阶段压缩空气储能工程项目的重点研发和示范方向。因此，从储气库方式来看，未来压缩空气储能电站储气方式将以人工硐库为主，其他储气方式为辅。

压缩空气储能电站经过5kW～1.5MW的实验性电站建设的验证，目前已投产项目基本达到了60MW～100MW等级，300MW等级项目也已在规划设计过程中。300MW等级电站储气库容积一般为10万m

在如此高压差变幅充放气过程中，硐库内空气的温度会随充放气过程产生明显变化，充气阶段温度升高，抽气阶段温度降低。硐库内空气温度的变化，除受空气热力学性质控制外，还受空气与硐壁之间的对流换热、以及硐壁各介质层的热传导过程的耦合控制。

在工程设计阶段，明晰硐库内空气压力和温度变化是发电机组设计的基础。另一方面，随着硐壁热传导过程的进行，密封层、混凝土衬砌以及围岩温度的变化，对其应力、应变和稳定性有重要影响。因此进行硐库热力学、硐壁对流换热以及硐壁热传导的耦合计算是一项重要的工作内容。

基于硐库内气体能量和质量守恒建立能量守恒方程和质量守恒方程，再根据热传导建立硐壁多层介质的热传导方程，以及各自的边界条件进行耦合求解，将是十分复杂的。

已有计算方法仅能计算标准循环条件下的充抽气过程，不能计算从大气压条件下开始充抽气的过程，而该过程对设计起到控制性作用，因此寻求更加适用的求解方法十分必要。

发明内容

本专利提出一种基于时间离散的耦合方法，可以方便地耦合硐库热力学以及硐壁热传导过程。

本发明所采用的技术手段如下所述：

一种地下硐库热力学与硐壁热传导耦合计算方法，包括以下步骤：

S1：确定硐库几何尺寸、运行参数，给出密封层、衬砌和围岩的热力学和力学基本计算参数，确定总计算时长t

S2：计算每个时步硐内空气充气和抽气速率、充气温度，以及硐库内空气密度；

S3：根据硐壁与空气之间的对流换热与硐壁热传导能量守恒，基于第i时步的硐库空气温度和硐壁温度分布，计算第i+1时步密封层的温度T

S4：假定密封层温度在第i+1时步中不变，计算第i+1时步硐库内空气的温度，计算硐壁热传导，得到硐壁温度分布；

S5：重复步骤S3～S4，直至大于总计算时长。

作为优选，在每一时步循环计算中，根据硐库空气与硐壁对流换热功率与硐壁表面热传导功率相等，求得密封层温度T

作为优选：

将拟计算的时域进行离散，选取时步Δt；

以硐库中心向外延伸设坐标轴r，并设坐标轴r的离散间隔为Δr；

假定时步内硐库密封层温度T

以此计算一时步结束时硐内气体温度T和围岩温度梯度

式中：

k为混凝土衬砌的导热系数，W·(m·K)

h

T

式中：

h

为围岩热传导功率，W；

进一步由式(2)计算得到下一时间间隔计算时需要的密封层温度T

作为优选，在计算硐内气体温度T基本解时，假设硐壁温度恒定，硐库内空气密度不变，充气速率或抽气速率为常数，以及空气为理想状态气体，进一步将下式(3)的能量守恒方程转化为常微分方程，

式中：

V为硐库体积，m

ρ为硐库内空气的密度，kg/m

c

T为硐库内空气的温度，K；

t为时间，s；

m

m

h

h为硐库内空气的焓，J；

Z为空气压缩系数；

R为空气常数，286.7J/(kg·K)；

u为硐库内空气的内能，J；

为硐库内空气与硐壁的换热速率，J/s；

并结合下式(4)的质量守恒方程和式(5)的对流换热方程得到硐内气体温度的解析解式(6)，

式中：

h

A

r

式(6)中：

其中：

ρ

ρ

T

m

m

V为硐库体积，m

T

作为优选，在计算围岩热传导基本解时，围岩热传导方程为下式(10)：

式中：

k

ρ

c

求解式(10)的初始条件和边界条件为：

T

T

T

式中：

r

f(r)为某计算开始时间点，围岩温度随计算点位置的函数，由上一时间段计算结果确定；

T

T

对式(10)两侧进行拉普拉斯变换，得：

式中：

s为拉普拉斯算子；

a

T

为T

整理式(14)得：

式(17)为虚宗变量零阶非齐次贝塞尔方程，其对应的齐次方程的通解为：

式中：

C

I

采用常数变易法，可以求得方式(17)的一个特解为：

由此得到式(17)的通解为：

将边界条件式(12)和式(13)进行拉普拉期变换后，代入式(20)，整理后得：

式中：

至此求得对流换热模式下，围岩中温度分布函数的普拉斯变换

围岩中温度分布函数的普拉斯变换

作为优选，通过下式(23)对式(22)进行求解：

式中，

作为优选，在对式(23)求解时，N取8～20之间的偶数。

作为优选，在对式(23)求解时，N取18，计算结果精度范围为10

作为优选，时步长度Δt取10s～600s。

与现有技术相比，本发明的有益技术效果在于，本耦合方法是在以密封层温度为边界条件的基础上将硐内空气热力学计算与围岩热传导分别进行，较以往根据热传导建立硐壁多层介质的热传导方程以及各自的边界条件进行耦合求解，本方法更为简便高效。

附图说明

图1为本方案考虑硐内热力学过程及硐壁热传导的储气库计算流程图。

图2为图1中每一循环步中的计算流程图。

图3为压缩空气储能电站技术原理图。

图4为隧道式地下人工硐库的纵剖图(a)和横剖图(b)。

图5为大罐式压缩空气储能密封硐库示意图。

图6为密封层、衬砌和围岩传热简化模型示意图。

图7为充气和抽气速率函数示意图。

图8为时间离散方法示意图。

具体实施方式

如图1所示，提出一种地下硐库热力学与硐壁热传导耦合计算方法，包括以下步骤：

S1：确定硐库几何尺寸、运行参数，给出密封层、衬砌和围岩的热力学和力学基本计算参数，确定总计算时长t

S2：计算每个时步硐内空气充气和抽气速率、充气温度，以及硐库内空气密度；

S3：根据硐壁与空气之间的对流换热与硐壁热传导能量守恒，基于第i时步的硐库空气温度和硐壁温度分布，计算第i+1时步密封层的温度T

S4：假定密封层温度在第i+1时步中不变，计算第i+1时步硐库内空气的温度，计算硐壁热传导，得到硐壁温度分布；

S5：重复步骤S3～S4，直至大于总计算时长。

如图2所示，在以上每一时步循环计算中，根据硐库空气与硐壁对流换热功率与硐壁表面热传导功率相等，求得密封层温度T

具体的，对于计算的物理模型：

由硐库的长度一般远大于其横向直径，可以将硐库的热传导过程作为轴对称问题考虑。

对于采用钢衬作为密封层的硐库，钢衬的厚度一般在10mm～20mm，厚度较薄。再考虑到其导热系数一般是衬砌和围岩的10倍以上，因此，可以不考虑钢衬的阻热，假定钢衬内任意一点的温度都相等。

再考虑到混凝土衬砌的导热系数、等压比热[2.94W/(m·K)，0.96kJ/(kg·K)]与围岩[2～5W/(m·K)，0.75～1.14kJ/(kg·K)]接近，将混凝土衬砌与围岩视为同一种传热介质，具有相同的热物性参数，本方案中将二者统称围岩。再考虑到硐内空气温度变化在衬砌和围岩中的传导深度有限，衬砌材料的参数对温度分布的影响大一些，所以混凝土衬砌与围岩材料参数均取混凝土衬砌的参数。

如图6所示，钢衬外侧半径记为r

以及对于计算的数学模型：

已有的均采用联立求解的方法求解以下能量守恒方程、质量守恒方程、对流换热方程、硐内空气与硐壁热传导边界条件、围岩热传导方程、广义气体状态方程。

其中，能量守恒方程为：

式中：

V为硐库体积，m

ρ为硐库内空气的密度，kg/m

c

T为硐库内空气的温度，K；

t为时间，s；

m

m

h

h为硐库内空气的焓，J；

Z为空气压缩系数；

R为空气常数，286.7J/(kg·K)；

u为硐库内空气的内能，J；

为硐库内空气与硐壁的换热速率，J/s。

式(3)中充气速率为正，抽气速率为负，具体定义见图7。图7还显示，在充气阶段和抽气阶段，其速率均为常数。

同时还有质量守恒方程：

对流换热方程：

式中：

h

A

r

硐内空气与硐壁热传导边界条件：

T

式中：

k为混凝土衬砌的导热系数，W·(m·K)

围岩热传导方程：

式中：

k

ρ

c

广义气体状态方程：

p＝ZρRT (25)

式中：

p为硐库内气体压力，kPa；

Z为空气压缩系数，无量纲；

R为空气常数，286.7J/(kg·K)。

在上述以往求解中，由于能量守恒方程的复杂性，其求解过程一般要假定一个周期内空气密度不变。当硐库工作压力变化范围较大时，这将带来较大的计算误差。已有解不能计算首次充气过程中硐库内空气的变化，而首次充气由于空气压力和密度变化大，其温度变化范围大，是工程设计的控制性因素。

因此对于上述已有计算方法中存在的缺陷，本方案通过：

将拟计算的时域进行离散，选取时步Δt，如图8所示，在每一个时间段内，假定硐库内空气密度不变，且密封层内表面温度不变；假定时步内硐库密封层温度T

且实际计算时，以硐库中心向外延伸设坐标轴r，将坐标轴r进行离散，并设坐标轴r的离散间隔为Δr；

以此为边界条件，计算一时步结束时硐内气体温度T和围岩温度梯度

式中：

k为混凝土衬砌的导热系数，W·(m·K)

h

T

式中：

h

为围岩热传导功率，W；

进一步由式(2)计算得到下一时间间隔计算时需要的密封层温度T

其中，在计算硐内气体温度T基本解时，假设硐壁温度恒定，硐库内空气密度不变，充气速率或抽气速率为常数，以及空气为理想状态气体(Z＝1)，进一步将下式(3)的能量守恒方程转化为常微分方程，

式中：

V为硐库体积，m

ρ为硐库内空气的密度，kg/m

c

T为硐库内空气的温度，K；

t为时间，s；

m

m

h

h为硐库内空气的焓，J；

Z为空气压缩系数；

R为空气常数，286.7J/(kg·K)；

u为硐库内空气的内能，J；

为硐库内空气与硐壁的换热速率，J/s；

并结合下式(4)的质量守恒方程和式(5)的对流换热方程得到硐内气体温度的解析解式(6)，

式中：

h

A

r

式(6)中：

其中：

ρ

ρ

T

m

m

V为硐库体积，m

T

以及，在计算围岩热传导基本解时，当密封层温度与硐内空气温度不等时，存在对流换热，围岩热传导方程为下式(10)：

式中：

k

ρ

c

求解式(10)的初始条件和边界条件为：

T

T

T

式中：

r

f(r)为某计算开始时间点，围岩温度随计算点位置的函数，由上一时间段计算结果确定；

T

T

对式(10)两侧进行拉普拉斯变换，得：

式中：

s为拉普拉斯算子；

a

T

为T

整理式(14)得：

式(17)为虚宗变量零阶非齐次贝塞尔方程，其对应的齐次方程的通解为：

式中：

C

I

采用常数变易法，可以求得方式(17)的一个特解为：

由此得到式(17)的通解为：

将边界条件式(12)和式(13)进行拉普拉期变换后，代入式(20)，整理后得：

式中：

至此求得对流换热模式下，围岩中温度分布函数的普拉斯变换

围岩中温度分布函数的普拉斯变换

考虑到温度函数计算时步比较短，通过下式(23)对式(22)进行数值计算求解：

式中，

在式(23)中，原则上，N取值越大，计算越精确，但是由于舍入误差的影响，在对式(23)求解时，N一般取8～20之间的偶数；当N取18时，计算结果精度范围达10

以及，时间步长越短，总计算耗时越长，且易出现数值计算困难；时间步长越长，计算误差越大，计算耗时越短。实际计算显示，时步长度Δt取10s～600s时，计算结果精度满足要求。