CN104504755B - 一种山体中分布式地下设施温度场仿真方法 - Google Patents

Abstract

本发明公开了一种在渗流效应影响下对含有分布式地下设施的山体进行温度场仿真的方法。利用高程信息提取的等高线数据建立山体和地下设施的三维几何模型，将渗流场等效成山体随机均匀分布的“毛细管”，将山体数据抽象成一棵具有层次结构的多叉树，利用计算机图形学里判断点在封闭图形的算法，精确计算出每根“毛细管”的高度，从而在山体几何模型中建立了等效渗流场的几何模型。再利用ANSYS生成的配置文件，通过程序化设计查找出构建的地下设施和“毛细管”的面的信息，对在渗流场影响下含有分布式地下设施的山体进行温度场的仿真。

Description

一种山体中分布式地下设施温度场仿真方法

技术领域

本发明属于自然地理、热物理学和信息处理交叉的领域，更具体地，涉及一种山体中分布式地下设施温度场仿真方法，用于在渗流效应影响下对含有分布式地下设施的山体进行温度场仿真。

背景技术

裂隙岩体是坝基、边坡、地下洞室等岩体工程中广泛遇到的一类复杂岩体，是由随机分布的裂隙和被裂隙切割的岩块组成的不连续介质，岩石与泥土中有空隙，富含有流动的水介质，称为渗流。

山体温度场的分布受地壳浅层温度分布和热状态的控制，也受岩体工程和外界条件(包括地下水渗流)的影响。理论和实验研究表明，热力的传递(即传热)一般有3种方式:传导、对流和辐射。工程岩体中温度场的分布是以传导和对流为主实现的，地下水本身的渗流运动通过热对流的方式进行了热能转移，从而影响了裂隙岩体的温度场分布；地下水在裂隙岩体中的渗流速度的大小直接控制了岩体温度的变化幅度，因此，渗流对温度的影响是非常明显的。

由于断裂带内存在各种裂隙及破碎的岩石，地下水流在裂隙中沿垂直方向流动，地下水流的温度低于岩体的温度，因此，岩体的温度传递给水流，使水流的温度在流动方向上逐渐升高，进而改变了岩体温度场的分布。通过水流与围岩的热对流交换，水流的温度近似于岩体的温度，在给定的范围内热量达到了平衡状态。

现有的地下物体热仿真方法，只考虑到地下设施在岩石土壤等固体传导下的热场分布，没有考虑到岩石裂隙及其中的流体介质对地下温度场分布的影响，而流体介质使得地下设施与山体的热场交换更快速，在表面表现得更明显，在大型通用有限元分析软件(ANSYS)中，也没有描述岩石裂隙的渗流率参数。

发明内容

本发明认为裂隙岩体是一种具有连续介质性质的物质，即假设裂隙岩体是由空隙性差而导水性强的裂隙系统和空隙性好而导水性弱的岩块所构成的，根据双重介质理论，本发明提出了一种将渗流场等效成山体随机均匀分布的“毛细管”的方法，对在渗流场影响下含有分布式地下设施的山体进行温度场的仿真。

为了实现上述目的，本发明提供了一种山体中分布式地下设施温度场仿真方法，所述方法包括如下步骤：

(1)山体及分布式地下设施的物理特征和温度分布仿真建模：对山体及分布式地下设施进行物理特征建模，根据物理特征模型在地下目标的温度/红外成像仿真平台(Simulation Platform Of Temperature/Infrared Image Formation，SPTIF)中直接对山体和分布式地下设施的几何结构进行物理特征建模，从而建立山体和分布式地下设施的几何模型；

(2)渗流场建模：将渗流场抽象为多根山体中随机均匀分布的细导管，加入到分布式地下设施温度仿真模型中；

(3)对包含分布式地下设施的山体做含有渗流场和不含渗流场的情况下进行温度场仿真。

在本发明的一个实施例中，所述步骤(2)具体包括如下子步骤：

(2.1)对山体数据的数据结构抽象，将整个山体抽象成一个有限节点组成的有层次关系的山体树结构；

(2.2)利用判断点是否在一个封闭多边形里的算法遍历上述生成的山体树结构，生成多根随机细管的底面坐标，得到每根细管的高度，生成毛细管。

在本发明的一个实施例中，所述步骤(2.1)具体包括：

将山体等高线数据分离成17个区域，每个区域在同一海拔高度上，不同区域之间在水平上相互独立，互不交错，在垂直方向具有层次结构，海拔高的区域在垂直方向上的投影包含于海拔低的封闭曲线内；

用多叉树的数据结构对山体进行进一步的抽象，山体树结构的构建方法是，把山体倒过来，每个区域当做一个节点，山底就是树的根节点，若在垂直方向上的投影中上一层对下一层有封闭区间的包含关系，那么下一层就是上一层的子节点。

在本发明的一个实施例中，所述步骤(2.2)具体包括：

随机产生一根“毛细管”的底面坐标，从山体树结构的最低层的节点开始，判断该坐标点是否在该层的所有点组成的封闭曲线内，若在此封闭曲线内，如果该层所在的节点是叶子节点，则“毛细管”的最大高度等于此节点的海拔高度，如果不是叶子节点，则对下一层子节点进行遍历；若不在此封闭曲线内，则对同一海拔的其它节点进行遍历，如果该点不在此海拔高度的所有节点的封闭曲线内，说明该“毛细管”的最大高度等于此节点的父节点的海拔高度；对山体树结构中的所有点遍历，找到所有点对应的最高海拔，即确定了所有“毛细管”对应的高度。

在本发明的一个实施例中，所述步骤(2.2)中判断点是否在一个封闭多边形里的算法具体包括：

(2.2.1)从P点作一条水平直线，从封闭多边形一点P₀开始，遍历完整个多边形所有的点，它前一个点表示为P₁，后一个点表示为P₂。

(2.2.2)求水平直线与多边形的所有交点，如果线段P₀P₂水平，如果p.y＝p₂.y，则把点P₂加入交点集合。如果线段P₀P₂水平，p₁.y＝p₀.y，则再次把点P₂加入交点集合。

(2.2.3)如果线段P₀P₂不水平，求出y＝p.y直线和线段P₀P₂的交点IP，如果IP和P₀重合，判断线段P₁P₀和线段P₀P₂是否在直线y＝p.y两侧，如果是，则把IP加入交点集合。如果IP不和P₀重合，则直接把IP加入交点集合。

(2.2.4)将交点集合点按横坐标大小排序；

(2.2.5)点在边界上判断为不在多边形内，若交点个数为奇数，判断点在多边形外，依次在交点集合中取两个点IP₁,IP₂，如果存在p.x＞＝p₁.x且p.x＜＝p₂.x，则点P在多边形里面，如果不存在，则点P在多边形外面。

在本发明的一个实施例中，所述步骤(3)具体包括：

(3.1)建立分布式地下设施热力学模型；

(3.2)含有分布式地下设施的山体温度场仿真。

在本发明的一个实施例中，所述步骤(3.1)具体包括：

将地下设施及其附近方形区域表示为Ω＝{x:0＜x_i＜l_i,i＝1,2,3}，方形区域沿三个维度x₁,x₂,x₃方向的边长为l₁,l₂,l₃，区域内任意一点的位置用x＝(x₁,x₂,x₃)表示，目标区域为Ω₁，背景区域为Ω\Ω₁，目标物体为一圆柱体，上表面位置为ρ₁，下表面位置为ρ₂，α_o和α_s分别为目标和背景区域的热扩散系数，k_o和k_s分别是目标和背景区域的热传导数数，观测持续时间为(0,t_e)，区域内任意一点的温度分布记为T(x,t),(x,t)∈Q_te＝Ω_x(0,t_e)；

区域Ω内任意一点的温度T(x,t)满足下面的偏微分方程：

在目标和周围土壤的交界面处，温度分布具有连续性，满足以下两个约束条件：

当时，n是空间沿x₁、x₂或者x₃方向的单位法向量；

其中，所述分布式地下设施热力学模型满足如下初始条件和边界条件，包括：

●初始条件：在起始观测时刻，假设地下岩土的温度分布已知，记为T(x,0)，

T(x,0)＝g(x),x∈Ω (5)

g(x)为温度分布函数，通过对起始观测时刻的岩土表面的温度分布和一些不同深度的温度样值进行插值处理得到；

●区域表面热量交换：区域起始观测时刻的温度分布已知后，需要确定其变化情况；给出一定的边界条件，利用导热微分方程求解出区域内部温度随空间和时间变化的分布状态，由此可以得出目标表面的温度分布，地物的边界条件表达为：

其中q_sun和q_sky分别为土壤吸收的太阳和天空辐射，q_conv为区域表面和空气之间通过热对流吸收的热量.经过数学变形，(8)式可以写成如下的线性形式

其中p和q(t)分别为天气条件和土壤热特性的函数，T_air是大气温度，T₀是土壤表面的温度分布，h_conv是土壤和大气之间的热对流系数。

●底面条件：假设在足够深的土壤处的温度分布是恒定不变的；

其中T_∞通过对一些位置的测量值进行插值处理得到；

●垂直边界条件：假设选取的岩土区域足够大，在方形区域除去上下表面后剩余四个面处土壤温度分布满足下面的边界条件

n是岩土区域Ω外表面除去上下表面后剩余的四个面上的向内或者向外的单位法向量；

根据方程(1)—(4)以及初始条件(5)和边界条件(8)(12)(13)构成了包含浅层地下设施区域温度分布的热物理模型，通过对该模型的求解，得到待预测区域表面的温度分布：

在本发明的一个实施例中，所述步骤(3.2)具体包括：在含有渗流场和不含渗流场的情况下对含有分布式地下设施的山体进行温度场仿真，在仿真过程中，将渗流场等效为多根“毛细管”，每根“毛细管”是一个六面的柱体，分布式地下设施也是由多段柱体模型组成，根据每根导管的每个面在集成了ANSYS的SPTIF中的编号以及导管的底面、顶面和侧面的信息，以及分布式地下设施的表面编号信息，从而对相应的表面赋温度，进行山体温度场分布仿真。

本发明方法通过对山体模型数据的数据结构抽象，精确计算山体随机均匀分布的“毛细管”的高度，构建出渗流场，再利用仿真软件的信息文件，通过自动检测算法，求解出“毛细管”和分布式地下设施的底面、顶面和侧面信息及其编号，对在渗流场影响下含有分布式地下设施的山体进行温度场的仿真。

附图说明

图1是本发明实施例中SPTIF平台主界面；

图2是本发明实施例中地下洞库几何结构平面图；

图3是本发明实施例中山体的三维立体图侧面视角；

图4是本发明实施例中山体等高线数据在ANSYS中的显示；

图5是本发明实施例中分布式地下设施的几何模型；

图6是本发明实施例中分布式地下设施的有限元模型；

图7是本发明实施例中山体的三维立体结构图；

图8是本发明实施例中山体的有限元模型；

图9是本发明实施例中山体树形结构；

图10是本发明实施例中判断点是否在一个封闭多边形里面示意图；

图11是本发明实施例中“毛细管”高度算法流程图；

图12是本发明实施例中添加“毛细管”后的山体剖面图；

图13是本发明实施例中地下物体及附近区域土壤的热模型示意图；

图14是本发明实施例中获取“毛细管”面信息和编号的算法流程图；

图15是本发明实施例中添加渗流现象地下设施三维建模；

图16是本发明实施例中无渗流现象山体热场仿真俯视图；

图17是本发明实施例中添加渗流现象山体热场仿真俯视图；

图18是本发明实施例中无渗流现象的山体热场仿真左视剖面图；

图19是本发明实施例中添加渗流现象的山体热场仿真左视剖面图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明利用以ANSYS软件作为计算地下设施温度场分布的内核需求，结合高级编程语言VC++和APDL语言二次开发地下目标温度/红外成像仿真平台(SPTIF)，建立分布式地下设施物理特征建模。

地下目标的温度/红外成像仿真平台(SPTIF)介绍：

地下目标的温度/红外成像仿真平台是针对本项目研究地下设施红外仿真图像的需要，以ANSYS软件作为计算地下设施温度场分布的内核需求，结合高级编程语言VC++和APDL语言二次开发的仿真平台。其核心思想是：借助ANSYS内核计算目标热量传递的过程，生成相应的温度分布图，然后根据物体的红外辐射原理计算其红外仿真图像。集成了ANSYS软件的内核以及VC的界面显示，使得热分析的过程更易掌握使用。并且可以在温度分布图的基础上，通过开发的相应的温度转红外算法将ANSYS解算的温度数据转换到红外数据，得到红外仿真图像。

地下目标温度/红外成像仿真平台如图1所示：

根据双重介质理论，将渗流场抽象为多根山体中随机均匀分布的很细的导管，进行山体分布式地下设施的温度场仿真，从而最大程度体现具有分布式地下设施的山体在其表面及内部的温度场分布情况，对指导地下设施的探测识别具有及其重要的意义，在现有的国内外文献中还没有看到与本发明相同的报道。

本发明提出了一种包含渗流效应的山体分布式地下设施温度场仿真方法，包括下述几个步骤：

(1)山体及分布式地下设施的物理特征和温度分布仿真建模：对山体及分布式地下设施进行物理特征建模，根据物理特征模型在SPTIF中直接对山体和分布式地下设施的几何结构进行物理特征建模，从而建立山体和分布式地下设施的几何模型；

(1.1)山体及分布式地下设施的物理特征建模

我们针对分布式地下设施特征(红外特征、多谱/高光谱特征、电磁波特征等)及其与环境的关系先期开展了理论分析，并针对对典型地下洞库几何结构利用MultiGen-Creator建立了目标物理模型，如图2所示。

此分布式地下设施是一个在山体中的地下洞库，利用Matlab软件绘制了地下洞库的自然场景-山体的三维立体图，其侧面视角如图3所示。

SPTIF几何模型的构建是以点、线、面和体来构成的，点(坐标)是构建几何模型的基础。几何体的创建是由关键点生成闭合的曲线，闭合的曲线生成平面，然后由闭合的曲面围成几何体。因此针对我们的山体三维立体图，我们需要山体表面的点坐标。我们通过Matlab中等高线的提取来获得山体某一高度的轮廓点坐标，在ANSYS中绘制等高线分布图如图4所示。

(1.2)山体及分布式地下设施温度分布仿真建模

根据物理特征模型，在SPTIF中直接建立山体和分布式地下设施的几何模型，其主要巷道是呈拱形状。山体的海拔高度约为160米，分布式地下设施高度为6米。

分布式地下设施几何模型如图5所示，其经过有限元网格划分后的三维模型如图6所示。

输入的山体等高线数据保存在17个不同的文件中，每个文件包含一组组成封闭曲线的所有点的信息，同一文件中的点所在的海拔高度相同，因此这些数据将整个山体分成了17个区域，这些区域在水平上相互独立，互不交错，在垂直方向具有层次结构，海拔高的区域在垂直方向上的投影包含于海拔低的封闭曲线内，这些数据结构特征对于渗流场建模的建立具有重要意义。

通过等高线数据将关键点连接成曲线，封闭曲线生成成平面，然后在侧面用蒙皮的构造形成构造结构体的封闭曲面进而实现复杂几何体的仿真实现，其三维立体结构图如图7所示，再经过有限元网格划分，三维模型如图8所示。

(2)渗流场建模，将渗流场抽象为多根山体中随机均匀分布的很细的导管，加入到分布式地下设施温度仿真模型中。

根据双重介质理论，将渗流场抽象为多根山体中随机均匀分布的很细的毛细导管，地下水在其中运动，运动的方向与多种因素有关，在这里我们仅仅把毛细管抽象成竖直方向，即水在垂直方向上运动。首先随机生成毛细导管的底面坐标，但是需要进一步确定每根毛细导管的高度，具体步骤如下：

(2.1)对山体数据的数据结构抽象，将整个山体抽象成一个有限节点组成的有层次关系的山体树结构；

将山体等高线数据分离成17个区域，每个区域在同一海拔高度上，不同区域之间在水平上相互独立，互不交错，在垂直方向具有层次结构，海拔高的区域在垂直方向上的投影包含于海拔低的封闭曲线内。按照海拔高度由高到低的顺序，我们将区域编号为1、4、5、6、(71、72)、(81、82)、9、10、(111、112、113)、(121、122、123)、13这17个编号，括号内的区域表示在同一海拔高度。

考虑到山体是分层构建，同一层的不同块相互独立，相邻两层的各块在空间上有一定的关系，整个山体是一个有限节点组成的一个有层次关系的集合，所以，选择用多叉树的数据结构对山体进行进一步的抽象，用来定位每根“毛细管”在空间上属于哪些块，从而确定随机分布的“毛细管”的高度信息。山体的树结构构建的大体思路是，把山体倒过来，每个区域当做一个节点(Node)，山底就是树的根节点，若在垂直方向上的投影中上一层对下一层有封闭区间的包含关系，那么下一层就是上一层的子节点，每个节点的信息如下：

构造山体的树形结构如图9所示，括号内表示该层区域对应的海拔高度。

(2.2)利用判断点是否在一个封闭多边形里的算法遍历上述生成的山体树结构，生成多根随机细管的底面坐标，得到每根细管的高度，生成毛细管。

要判断垂直方向的投影的包含关系要用到计算机图形学里的判断点是否落在封闭图形里的知识，本发明设计的算法原理如下：

这里用到了水平/垂直交叉点判别法(适用于任意封闭多边形，包括凹多边形和凸多边形)。如图10所示，从P点作水平直线，求解出直线与多边形的所有交点，在x轴方向，从左到右依次取出两个交点IP₁、IP₂，如果P在多边形内部，则必然存在点P在IP₁、IP₂两点的中间(若重合也算作中间)。所以，我们可以顺序考虑多边形的每条边，求出交点的总个数。

在求解直线和多边形交点的过程中，存在一些特殊情况，本算法约定如下，可以正确判断点是否在封闭多边形里面。对于多边形上一点P₀，考虑边P₀P₂和P₀的前一个点P₁：

a、如果只有P₀P₂水平，当P在P₀P₂所在直线上时，把P₂加入交点集合；

b、如果线段P₁P₀，P₀P₂都水平，当P₁也在P₀P₂所在直线上时，则两次把P₂加入交点集合。

c、如果y＝p.y直线与其交点是P₀，如果P₁P₀，P₀P₂都在y＝p.y的同一侧，则两次把P₀加入交点集合。

判断一点P是否在封闭多边形里面，p.x，p.y分别表示点P的横纵坐标，其具体算法流程如下：

(2.2.1)从P点作一条水平直线，从封闭多边形一点P₀开始，遍历完整个多边形所有的点，它前一个点表示为P₁，后一个点表示为P₂。

(2.2.2)求水平直线与多边形的所有交点，考虑到特殊情况。如果线段P₀P₂水平，如果p.y＝p₂.y，则把点P₂加入交点集合。如果线段P₀P₂水平，p₁.y＝p₀.y，则再次把点P₂加入交点集合。

(2.2.3)如果线段P₀P₂不水平，求出y＝p.y直线和线段P₀P₂的交点IP，如果IP和P₀重合，判断线段P₁P₀和线段P₀P₂是否在直线y＝p.y两侧，如果是，则把IP加入交点集合。如果IP不和P₀重合，则直接把IP加入交点集合。

(2.2.4)将交点集合点按横坐标大小排序

(2.2.5)点在边界上判断为不在多边形内，若交点个数为奇数，判断点在多边形外，依次在交点集合中取两个点IP₁,IP₂，如果存在p.x＞＝p₁.x且p.x＜＝p₂.x，则点P在多边形里面，如果不存在，则点P在多边形外面。

(2.3)生成多根随机细管的底面坐标，得到每根细管的高度，生成毛细管：

利用构建的“山体树”和判断点是否在一个封闭多边形里的算法遍历，其算法流程图如图11所示。我们按照此算法对山体遍历，找到对应某点的最高海拔，也就是确定了该“毛细管”的高度。

当随机产生一根“毛细管”的底面坐标，从最低层的节点开始(图9中的结点1)，判断该坐标点是否在该层的所有点组成的封闭曲线内，若在此封闭曲线内，如果该层所在的节点是叶子节点，则“毛细管”的最大高度等于此节点的海拔高度，如果不是叶子节点，则对下一层子节点进行遍历；若不在此封闭曲线内，则对同一海拔的其它节点进行遍历；如果该点不在此海拔高度的所有节点的封闭曲线内，说明该“毛细管”的最大高度等于此节点的父节点的海拔高度，；对山体树结构中的所有点遍历，找到所有点对应的最高海拔，即确定了所有“毛细管”对应的高度。

利用SPTIF进行仿真得到喻家山体，在山体中添加“毛细管”后的山体剖面图如图12所示。

(3)对包含分布式地下设施的山体做含有渗流场和不含渗流场的情况下进行温度场仿真。

(3.1)建立分布式地下设施热力学模型

为了建立分布式地下设施的热模型，我们需要先做以下假设：

1)埋藏目标的区域介质(岩石和土壤)同质，且各向同性；

2)在一段时间内，土壤的含水量等物理条件是不变的；

3)物体完全埋藏于土壤岩石里。

理论分析中我们用一个地下设施及其附近区域的示意图作说明，如13所示，整个方形区域为Ω＝{x:0＜x_i＜l_i,i＝1,2,3}，方形区域沿x₁,x₂,x₃方向的边长为l₁,l₂,l₃区域内任意一点的位置用x＝(x₁,x₂,x₃)表示，目标区域为Ω₁，背景区域为Ω\Ω₁，目标物体为一圆柱体，上表面位置为ρ₁，下表面位置为ρ₂，α_o和α_s分别为目标和背景区域的热扩散系数，k_o和k_s分别是目标和背景区域的热传导数数。观测持续时间为(0,t_e)，区域内任意一点的温度分布记为T(x,t),(x,t)∈Q_te＝Ω_x(0,t_e)。

区域Ω内任意一点的温度T(x,t)满足下面的偏微分方程：

在目标和周围土壤的交界面处，温度分布具有连续性，可以得到以下两个约束条件:

当时，n是空间沿x₁、x₂或者x₃方向的单位法向量。

为了求出方程组(1)～(4)的解，我们还需要一些初始条件和边界条件，下面我们将详细描述如何在该模型中确定初始条件和边界条件。

●初始条件：和解决其它动态问题类似，我们需要知道在起始观测时刻的区域温度分布。在起始观测时刻，我们假设地下岩土的温度分布已知，记为T(x,0)

T(x,0)＝g(x),x∈Ω (5)

g(x)为温度分布函数；在实际处理问题的过程中，起始观测时刻的地下岩土温度分布是无法直接获得的，我们可以对起始观测时刻的岩土表面的温度分布和一些不同深度的温度样值进行插值处理得到。

●区域表面热量交换：区域起始观测时刻的温度分布已知后，我们还需要确定其变化情况。引起区域温度变化的主要因素是区域表面和外界环境的热量交换，其是通过热辐射、对流和热传导三种形式进行的。区域表面和外界环境之间的热量交换主要是通过热辐射和热对流进行的，热传导主要在区域表面和区域内部之间进行。区域表面从外部环境接收到的辐射主要来自太阳和天空中的大气，同时区域表面也在向外辐射能量。

a)太阳的辐射(q_sun)

区域接收到的太阳辐射能量主要包括太阳直射、太阳散射和地面反射三部分。区域表面接收到的太阳直接辐射能量与区域表面的发射率、日地之间的距离、太阳常数、大气质量、大气透明度、区域表面对太阳直接辐射的接收角系数有关。

b)天空大气的辐射(q_sky)

天空中大气的辐射也对区域温度的变化产生影响。大气辐射主要是一种长波辐射，大气在吸收了太阳热量和地球热量后，具有了一定的温度，所以也在时时刻刻的向外辐射能量。

c)热对流(q_conv)

通过区域表面和大气之间的对流换热过程进入区域的能量和区域温度、大气温度以及对流换热系数h(流体种类、流体流动方式、流体状态、固体表面几何形状和温差)有关。

给出一定的边界条件，利用导热微分方程可以求解出区域内部温度随空间和时间变化的分布状态，由此可以得出目标表面的温度分布。导热方程是以能量守恒定律和傅氏定律为基础。设密度为ρ；地物的比热容为c；t表示时间；k表示地物的导热系数；Θ_v表示发热功率。

边界条件描述了导热物体和其所处的周围环境之间在两者的边界处所进行的热交换。确定地物的边界条件对我们分析非稳态导热过程具有重要的意义。

式中:为地物和周围环境边界面上的向外方向的法线。

在该区域模型中上式可以写成如下形式：

其中q_sun和q_sky分别为土壤吸收的太阳和天空辐射，q_conv为区域表面和空气之间通过热对流吸收的热量.经过数学变形，(8)式可以写成如下的线性形式

其中p和q(t)分别为天气条件和土壤热特性的函数，T_air是大气温度，T₀是土壤表面的温度分布，h_conv是土壤和大气之间的热对流系数。

●底面条件：我们假设在足够深的土壤处的温度分布是恒定不变的。

其中T_∞可以通过对一些位置的测量值进行插值处理得到，对于大多数的岩土和气象条件，岩土深度一般取0.5m就可以满足该深度的土壤温度基本恒定。

●垂直边界条件：我们假设选取的岩土区域足够大，在方形区域除去和后剩余四个面处土壤温度分布满足下面的边界条件

n是岩土区域Ω外表面除去和剩余的四个面上的向内或者向外的单位法向量。这可以理解为地下设施对周围区域温度分布的影响范围是有限的。

方程(1)—(4)以及初始条件(5)和边界条件(8)(12)(13)构成了我们建立的求解包含浅层地下设施区域温度分布的热物理模型，通过对该模型的求解，可以准确预测区域表面的温度分布。

(3.2)含有分布式地下设施的山体温度场仿真

在含有渗流场和不含渗流场的情况下对含有分布式地下设施的山体进行温度场仿真，在仿真过程中，将渗流场等效为多根“毛细管”，每根“毛细管”是一个六面的长方体，分布式地下设施也是由多段长方体模型组成，我们需要精确知道每根导管的每个面在集成了ANSYS的SPTIF中的编号以及底面、顶面和侧面的信息，同时也需要知道分布式地下设施的表面编号信息，从而对相应的表面赋温度，进行山体温度场分布仿真。由于毛细管数量巨大，本发明设计了一种程序化自动检测长方体型的“毛细管”和分布式地下设施的底面、顶面和侧面信息及其编号的方法，具体步骤如下：

1)对SPTIF中含有渗流场和分布式地下设施山体模型的点、线、面、体信息进行数据处理

山体及分布式地下设施温度分布仿真建模时，输入的山体等高线数据保存的是一块块封闭曲线的点的坐标信息，通过等高线数据进行关键点连接成曲线，封闭曲线生成成平面，然后在侧面用蒙皮的构造形成构造结构体的封闭曲面进而实现复杂几何体。在集成了ANSYS的SPTIF配置文件中，分别存在模型的点、线、面、体的信息文件，其中，点的信息文件KLIST.txt包含了所有点的坐标、编号和其它信息，线的信息文件LLIST.txt包含了线的编号、两端点的编号、线的长度和其它信息，面的信息文件ALIST.txt包含了面的编号，围成面的线段编号等信息，体的信息文件VLIST.txt包含了体的编号、形成体的面编号等信息。

我们对这些信息文件进行处理，读取并存入为数据结构，其具体数据结构如下：

点的结构体如下表：

NO	点的编号
		x	点的横坐标
y	点的纵坐标
		z	点的海拔高度
flag	标记位，是否为感兴趣的点

线的结构体如下表：

NO	线的编号
		point[2]	线的两个端点结构体
flag	标记位，是否为感兴趣的线

面的结构体如下表：

NO	面的编号
		line[4]	围成封闭面的四条线结构体
flag	标记位，是否为感兴趣的面
		keypoint	面中心点

isbeyondtunnel

标记位，是否为地下设施上方的面

体的结构体信息如下表

NO	体的编号
		areas[6]	体的六个面结构体
bottom	体的底面
		up	体的上表面
sides[4]	体的四个侧面
		flag	标记位，是否为感兴趣的体
height	体的高度

在获取了点、线、面、体的数据结构后，具体检测长方体型的“毛细管”和分布式地下设施的底面、顶面和侧面信息及其编号的算法流程图如图14所示。

2)利用“毛细管”和分布式地下设施的面的信息，求解温度场分布

利用“毛细管”和分布式地下设施的面的信息，建立含有渗流现象的地下设施三维建模，如图15所示。

设定边界条件，加载热载荷。具体的物性参数如下表所示：

环境温度	35℃
		地层温度	25℃
目标生热率	5×101(W/m3)
		介质传导系数	3.2(W/m℃)
表面换热系数	2(W/m℃)

地库所在地区的经纬度为：北纬30.6°，东经114.1°，海拔180米。取7月21日上午9时，经过太阳辐射特性公式的计算，得出地球表面的法向太阳直射辐射照度为82.4W/m²，正南垂直立面上的太阳辐射照度为11.4W/m²。

我们综合考虑山体表面的不平滑度以及山体相对太阳的偏斜角，在向阳面取太阳辐射的照度为10W/m²，背阴面取太阳辐射的照度为8W/m²。即太阳辐射转化成热流密度分别为向阳面10W/m²，背阴面8W/m²。实际应用中根据目标所处位置实地考察测量为准。

分别对分布式地下设施在没有渗流场和有渗流场的情况下进行温度场分布仿真，得到如下图结果。未加渗流的山体热场仿真俯视图如图16所示，有渗流的山体热场仿真俯视图如图17所示。未加渗流的山体热场仿真左视剖面图如图18所示，添加渗流的山体热场仿真左视剖面图如图19所示。

上述实验说明：渗流使地下设施的热能向山表面传导的更快、表现的更明显。

实验结论：从实验中可以看出，地下水本身的渗流运动通过热对流的方式进行了热能转移，这些因素影响了裂隙岩体的温度场分布，，地下水在裂隙岩体中的渗流速度的大小直接控制了岩体温度的变化幅度，因此，渗流对温度的影响是非常明显的。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种山体中分布式地下设施温度场仿真方法，其特征在于，所述方法包括如下步骤：

(1)山体及分布式地下设施的物理特征和温度分布仿真建模：对山体及分布式地下设施进行物理特征建模，根据物理特征模型在SPTIF中直接对山体和分布式地下设施的几何结构进行物理特征建模，从而建立山体和分布式地下设施的几何模型；(2)渗流场建模：将渗流场抽象为多根山体中随机均匀分布的细导管，加入到分布式地下设施温度仿真模型中；具体包括：

(2.1)对山体数据的数据结构抽象，将整个山体抽象成一个有限节点组成的有层次关系的山体树结构；

(2.2)利用判断点是否在一个封闭多边形里的算法遍历上述生成的山体树结构，生成多根随机细管的底面坐标，得到每根细管的高度，生成毛细管；具体为：随机产生一根“毛细管”的底面坐标，从山体树结构的最低层的节点开始，判断该坐标点是否在该层的所有点组成的封闭曲线内，若在此封闭曲线内，如果该层所在的节点是叶子节点，则“毛细管”的最大高度等于此节点的海拔高度，如果不是叶子节点，则对下一层子节点进行遍历；若不在此封闭曲线内，则对同一海拔的其它节点进行遍历，如果该点不在此海拔高度的所有节点的封闭曲线内，说明该“毛细管”的最大高度等于此节点的父节点的海拔高度；对山体树结构中的所有点遍历，找到所有点对应的最高海拔，即确定了所有“毛细管”对应的高度；

(3)对包含分布式地下设施的山体做含有渗流场和不含渗流场的情况下进行温度场仿真。

2.如权利要求1所述的方法，其特征在于，所述步骤(2.1)具体包括：

将山体等高线数据分离成17个区域，每个区域在同一海拔高度上，不同区域之间在水平上相互独立，互不交错，在垂直方向具有层次结构，海拔高的区域在垂直方向上的投影包含于海拔低的封闭曲线内；

用多叉树的数据结构对山体进行进一步的抽象，山体树结构的构建方法是，把山体倒过来，每个区域当做一个节点，山底就是树的根节点，若在垂直方向上的投影中上一层对下一层有封闭区间的包含关系，那么下一层就是上一层的子节点。

3.如权利要求1所述的方法，其特征在于，所述步骤(2.2)中判断点是否在一个封闭多边形里的算法具体包括：

(2.2.1)从P点作一条水平直线，从封闭多边形一点P₀开始，遍历完整个多边形所有的点，它前一个点表示为P₁，后一个点表示为P₂；

(2.2.2)求水平直线与多边形的所有交点，如果线段P₀P₂水平，如果p.y＝p₂.y，则把点P₂加入交点集合；如果线段P₀P₂水平，p₁.y＝p₀.y，则再次把点P₂加入交点集合；

(2.2.3)如果线段P₀P₂不水平，求出y＝p.y直线和线段P₀P₂的交点IP，如果IP和P₀重合，判断线段P₁P₀和线段P₀P₂是否在直线y＝p.y两侧，如果是，则把IP加入交点集合；如果IP不和P₀重合，则直接把IP加入交点集合；

(2.2.4)将交点集合点按横坐标大小排序；

(2.2.5)点在边界上判断为不在多边形内，若交点个数为奇数，判断点在多边形外，依次在交点集合中取两个点IP₁,IP₂，如果存在p.x＞＝p₁.x且p.x＜＝p₂.x，则点P在多边形里面，如果不存在，则点P在多边形外面。

4.如权利要求1、2或3所述的方法，其特征在于，所述步骤(3)具体包括：

(3.1)建立分布式地下设施热力学模型；

(3.2)含有分布式地下设施的山体温度场仿真。

5.如权利要求4所述的方法，其特征在于，所述步骤(3.1)具体包括：

将地下设施及其附近方形区域表示为Ω＝{x:0＜x_i＜l_i,i＝1,2,3}，方形区域沿三个维度x₁,x₂,x₃方向的边长为l₁,l₂,l₃，区域内任意一点的位置用x＝(x₁,x₂,x₃)表示，目标区域为Ω₁，背景区域为Ω\Ω₁，目标物体为一圆柱体，上表面位置为ρ₁，下表面位置为ρ₂，α_o和α_s分别为目标和背景区域的热扩散系数，k_o和k_s分别是目标和背景区域的热传导数数，观测持续时间为(0,t_e)，区域内任意一点的温度分布记为T(x,t),(x,t)∈Q_te＝Ω_x(0,t_e)；

区域Ω内任意一点的温度T(x,t)满足下面的偏微分方程：

$\frac{\∂T(x,t)}{\∂t}={\mathrm{\α}}_o\underset{j=1}{\overset{3}{\Σ}}\frac{{\∂}^{2}T(x,t)}{\∂x_j^2},(x,t)\∈{\mathrm{\Ω}}_1\×(0,t_e)---\left(1\right)$

$\frac{\∂T(x,t)}{\∂t}={\mathrm{\α}}_s\underset{j=1}{\overset{3}{\Σ}}\frac{{\∂}^{2}T(x,t)}{\∂x_j^2},(x,t)\∈(\mathrm{\Ω}\backslash {\mathrm{\Ω}}_1)\×(0,t_e)---\left(2\right)$

在目标和周围土壤的交界面处，温度分布具有连续性，满足以下两个约束条件：

$\underset{y\∈{\mathrm{\Ω}}_1,y\→x}{\lim }T(y,t)=\underset{y\∈\mathrm{\Ω}\backslash {\mathrm{\Ω}}_1,y\→x}{\lim }T(y,t)---\left(3\right)$

$k_o\frac{\∂T{|}_{{\mathrm{\Ω}}_1}}{{\∂}_n}(x,t)=k_s\frac{\∂T{|}_{\mathrm{\Ω}\backslash {\mathrm{\Ω}}_1}}{{\∂}_n}(x,t)---\left(4\right)$

当时，n是空间沿x₁、x₂或者x₃方向的单位法向量；

其中，所述分布式地下设施热力学模型满足如下初始条件和边界条件，包括：

●初始条件：在起始观测时刻，假设地下岩土的温度分布已知，记为T(x,0)，

T(x,0)＝g(x),x∈Ω (5)

g(x)为温度分布函数，通过对起始观测时刻的岩土表面的温度分布和一些不同深度的温度样值进行插值处理得到；

●区域表面热量交换：区域起始观测时刻的温度分布已知后，需要确定其变化情况；给出一定的边界条件，利用导热微分方程求解出区域内部温度随空间和时间变化的分布状态，由此可以得出目标表面的温度分布，地物的边界条件表达为：

$-k_s\frac{\∂T}{\∂x_3}(x,t)=q_{sun}\left(t\right)+q_{sky}\left(t\right)+q_{conv}(x,t),(x,t)\∈{\mathrm{\Γ}}_3^1\×(0,t_e)---\left(8\right)$

其中q_sun和q_sky分别为土壤吸收的太阳和天空辐射，q_conv为区域表面和空气之间通过热对流吸收的热量.经过数学变形，(8)式可以写成如下的线性形式

$-{\mathrm{\α}}_s\frac{\∂T(x,t)}{\∂x_3}(x,t)+pT(x,t)\∈{\mathrm{\Γ}}_3^1\×(0,t_e)---\left(9\right)$

其中p和q(t)分别为天气条件和土壤热特性的函数，T_air是大气温度，T₀是土壤表面的温度分布，h_conv是土壤和大气之间的热对流系数；

$p=\frac{{\mathrm{\α}}_s}{k_s}{h}_{conv}---\left(10\right)$

$q\left(t\right)=\frac{{\mathrm{\α}}_s}{k_s}\[q_{sun}\left(t\right)+q_{sky}\left(t\right)+h_{conv}{T}_{air}\left(t\right)\]---\left(11\right)$

●底面条件：假设在足够深的土壤处的温度分布是恒定不变的；

$T(x,t)=T_{\mathrm{\∞}},(x,t)\∈{\mathrm{\Γ}}_3^2\×(0,t_e)---\left(12\right)$

其中T_∞通过对一些位置的测量值进行插值处理得到；

●垂直边界条件：假设选取的岩土区域足够大，在方形区域除去上下表面后剩余四个面处土壤温度分布满足下面的边界条件

$\frac{\∂T}{\∂n}(x,t)=0---\left(13\right)$

n是岩土区域Ω外表面除去上下表面后剩余的四个面上的向内或者向外的单位法向量；

根据方程(1)—(4)以及初始条件(5)和边界条件(8)(12)(13)构成了包含浅层地下设施区域温度分布的热物理模型，通过对该模型的求解，得到待预测区域表面的温度分布：

$\left\{\begin{array}{c}\frac{\∂T(x,t)}{\∂t}={\mathrm{\α}}_o\underset{i=1}{\overset{3}{\Σ}}\frac{{\∂}^{2}T(x,t)}{\∂x_i^2},(x,t)\∈{\mathrm{\Ω}}_1\×(0,t_e)\\ \frac{\∂T(x,t)}{\∂t}={\mathrm{\α}}_s\underset{i=1}{\overset{3}{\Σ}}\frac{{\∂}^{2}T(x,t)}{\∂x_i^2},(x,t)\∈(\mathrm{\Ω}\backslash {\mathrm{\Ω}}_1)\×(0,t_e)\\ T(x,0)=g\left(x\right),x\∈\mathrm{\Ω}\\ -k_s\frac{\∂T}{\∂x_3}(x,t)=q_{sun}\left(t\right)+q_{sky}\left(t\right)+q_{conv}(x,t),(x,t)\∈{\mathrm{\Γ}}_3^1\×(0,t_e)\\ T(x,t)=T_{\mathrm{\∞}},(x,t)\∈{\mathrm{\Γ}}_3^2\×(0,t_e)\\ \frac{\∂T}{\∂n}(x,t)=0\end{array}\right..$

6.如权利要求5所述的方法，其特征在于，所述步骤(3.2)具体包括：在含有渗流场和不含渗流场的情况下对含有分布式地下设施的山体进行温度场仿真，在仿真过程中，将渗流场等效为多根“毛细管”，每根“毛细管”是一个六面的柱体，分布式地下设施也是由多段柱体模型组成，根据每根导管的每个面在集成了ANSYS的SPTIF中的编号以及导管的底面、顶面和侧面的信息，以及分布式地下设施的表面编号信息，从而对相应的表面赋温度，进行山体温度场分布仿真。