CN109299502B - 一种连续-非连续介质热传导的二维数值模拟方法及系统 - Google Patents

Abstract

本发明提供了一种连续‐非连续介质热传导的二维数值模拟方法及系统，首先将连续‑非连续介质进行网格划分，在这些实体单元的网格的公共边上插入节理单元或者界面单元，或者不插入，标记节理单元、界面单元或实体单元的公共边；一个实体单元内的温度分布根据实体单元节点的温度进行插值来表示，实体单元和相邻实体单元通过节理单元或者界面单元或实体单元的公共边发生热量传递；根据实体单元内的热传导和实体单元间的热量交换计算节点的总热流量；根据有限差分积分策略，以节点的总热量不断更新节点的温度，最终完成连续‑非连续介质的热传导。本发明的有益效果是：本发明提供的技术方案操作简单，适用性和实用性强。

Description

一种连续-非连续介质热传导的二维数值模拟方法及系统

技术领域

本发明涉及能源环境、岩土工程领域，尤其涉及一种连续-非连续介质热传导的二维数值模拟方法及系统。

背景技术

连续-非连续介质热传导问题涉及地热开采、石油天然气开采、核废料处置等工程应用问题。随着计算机硬件的发展及计算仿真技术的发展，采用数值模拟方法逐渐成为传热分析计算的有力工具。目前，求解传热问题的主要计算方法是有限元和有限差分方法两种。这两种方法的基本原理是，将连续体划分网格，将整个连续体温度的求解转变为单元节点温度的求解。对于连续体的传热分析而言，这两种方法均有很好的适用性，可以处理各种边界情形和不同条件的传热问题，但对于连续-非连续介质，因为这两种方法没有考虑到裂缝对传热的阻隔作用，尤其是当介质中存在裂缝的动态扩展时，裂缝对传热的阻隔作用较大，更应该考虑到裂缝对传热的阻隔作用，所以有限元和有限差分这两种方法难以处理传热分析问题，而在工程和自然界中，有很多介质(如岩体)是连续-非连续介质，因此，亟需一种方法能够处理连续-非连续介质的传热分析问题。

发明内容

为了解决上述问题，本发明提供了一种连续-非连续介质热传导的二维数值模拟方法及系统，一种连续-非连续介质热传导的二维数值模拟方法及系统，主要包括以下步骤：

S101：对连续-非连续介质区域进行网格划分；网格划分所形成的网格为实体单元，网格的各个顶点作为节点；

S102：判断是否在网格间的公共边上插入了节理单元或者界面单元？若是，则到步骤S103，若否，则到步骤S104；公共边分为节理和裂隙非连续面处公共边和连续面处公共边；

S103：将插入在包括节理和裂隙非连续面处公共边的节理单元或者界面单元标记为断裂单元，将插入在连续面处公共边的节理单元或者界面单元标记为非断裂单元，然后到步骤S105；

S104：将位于包括节理和裂隙非连续面处的公共边标记为断裂边，将连续面处的公共边标记为非断裂边，然后到步骤S105；

S105：通过所述断裂单元或者非断裂单元、断裂边或者非断裂边，实体单元内发生热传导，相邻实体单元间发生热量传递；

S106：根据实体单元内的热传导和相邻实体单元间的热量传递，分别计算各个节点的总热流量；

S107：根据有限差分积分策略，利用计算得到的所述各个节点的总热流量分别不断更新对应节点的温度，最终完成连续-非连续介质区域的热传导。

进一步地，在步骤S106中，计算任意一节点的总热流量时，断裂单元和断裂边的换热系数需要乘以一个折减系数；根据相邻实体单元间的温度差和断裂单元或断裂边的乘以折减系数的换热系数，计算得到实体单元间经过所述任意一节点传递的热流量；根据实体单元的所有节点的温度，求得实体单元内的温度梯度，然后求得实体单元内流入所述任意一节点的热流量；将上述经过所述任意一节点传递的热流量和流入所述任意一节点的热流量相加得到所述任意一节点的总热流量。

进一步地，在步骤S107中，所述任意一节点的温度可按照以下公式不断进行更新：

其中，T^t+Δt是所述任意一节点在下一个时间步的温度，T^t是所述任意一节点在当前时间步t的温度，Δt是时间步长，C_p是比热容，M是所述任意一节点的质量。

一种连续-非连续介质热传导的二维数值模拟系统，其特征在于：包括网格划分模块、判断模块、插入标记模块、非插入标记模块、计算模块和更新模块；

网格划分模块，用于对连续-非连续介质区域进行网格划分；网格划分所形成的网格为实体单元，网格的各个顶点作为节点；

判断模块，用于判断是否在网格间的公共边上插入了节理单元或者界面单元？

插入标记模块，用于将插入在包括节理和裂隙非连续面处公共边的节理单元或者界面单元标记为断裂单元，将插入在连续面处公共边的节理单元或者界面单元标记为非断裂单元；

非插入标记模块，用于将位于包括节理和裂隙非连续面处的公共边标记为断裂边，将连续面处的公共边标记为非断裂边；

计算模块，用于分别计算各个节点的总热流量；计算任意一节点的总热流量时，所述断裂单元和断裂边的换热系数需要乘以一个折减系数；根据相邻实体单元间的温度差和所述断裂单元或断裂边的乘以折减系数的换热系数，计算得到实体单元间经过所述任意一节点传递的热流量；根据实体单元的所有节点的温度，求得实体单元内的温度梯度，然后求得实体单元内流入所述任意一节点的热流量；将上述经过所述任意一节点传递的热流量和流入所述任意一节点的热流量相加得到所述任意一节点的总热流量；

更新模块，用于计算得到的所述各个节点的总热流量分别不断更新对应节点的温度，最终完成连续-非连续介质区域的热传导。

进一步地，在更新模块中，所述任意一节点的温度可按照以下公式不断进行更新：

其中，T^t+Δt是所述任意一节点在下一个时间步的温度，T^t是所述任意一节点在当前时间步t的温度，Δt是时间步长，C_p是比热容，M是所述任意一节点的质量。

本发明提供的技术方案带来的有益效果是：本发明将连续-非连续介质进行网格区域化为实体单元，使计算过程更为容易理解、也更简单；通过乘以折减系数，考虑到裂缝对传热的阻碍作用和在裂缝的动态扩展过程中裂缝对热传导的影响，提高了计算结果的准确性，因此本发明提供的技术方案操作简单，适用性和实用性强。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种连续-非连续介质热传导的二维数值模拟方法及系统的流程图；

图2是本发明实施例中连续-非连续介质网格离散示意图；

图3是本发明实施例中插入节理单元或界面单元示意图；

图4是本发明实施例中标记节理单元、界面单元或实体单元公共边示意图；

图5是本发明实施例中计算示连续介质热传导示意图；

图6是本发明实施例中计算示连续介质热传导三角形单元通过节理单元、界面单元或实体单元公共边的热量交换示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种连续-非连续介质热传导的二维数值模拟方法及系统。本发明首先将没有作标记的节理单元、界面单元或实体单元的公共边的换热系数取一个较大的值，使得网格离散元化不影响连续体的传热性质；然后根据实体单元内节点的温度，求得实体单元内的温度梯度，然后求得实体单元本身流入各个节点的热流量；根据节理单元或者界面单元连接相邻实体单元间的温度差，以及节理单元或者界面单元的换热系数，计算经由节理单元或者界面单元流入节点的热流量；两个热流量相加即获得节点的总热流量，然后根据有限差分积分策略更新节点的温度；如此循环即可获得连续-非连续介质内温度的演化。

请参考图1，图1是本发明实施例中一种连续-非连续介质热传导的二维数值模拟方法及系统的流程图，具体包括如下步骤：

S101：对连续-非连续介质区域进行网格划分；网格划分所形成的网格为实体单元，网格的各个顶点作为节点；以连续-非连续介质的包括节理和裂隙的非连续面作为约束条件，可以采用现有商用有限元软件、前后处理软件以及一些开源网格程序对连续-非连续介质进行网格划分，将所述连续-非连续介质划分为包括三角形单元、四边形单元、voronoi单元和其他任意多变形中一种或多种实体单元的网格；所述实体单元由节点构成；比如三角形单元由3个节点构成；如图2为针对连续-非连续介质划分的网格；

S102：判断是否在网格间的公共边上插入了节理单元或者界面单元？若是，则到步骤S103，若否，则到步骤S104；公共边分为节理和裂隙非连续面处公共边和连续面处公共边；

S103：将插入在包括节理和裂隙非连续面处公共边的节理单元或者界面单元标记为断裂单元，将插入在节理和裂隙连续面处公共边的节理单元或者界面单元标记为非断裂单元，然后到步骤S105；如果在模拟的过程中有新的裂缝生成，即有新的节理单元或者界面单元被插入和标记为断裂单元；插入了节理单元或者界面单元后，网格的公共边进行分离成两条边，最终形成如图3所示的网格；

S104：将位于包括节理和裂隙非连续面处的公共边标记为断裂边，将连续面处的公共边标记为非断裂边，然后到步骤S105；如果在模拟的过程中有新的裂缝生成，即有新的实体单元的公共边被标记为断裂边；

标记位于节理、裂隙等非连续面的节理单元、界面单元或实体单元公共边的示意图如图4所示；

S105：通过所述断裂单元或者非断裂单元、断裂边或者非断裂边，实体单元内发生热传导，相邻实体单元间发生热量传递；一个实体单元内的温度分布是通过对实体单元的温度进行插值来表示的；

S106：根据实体单元内的热传导和相邻实体单元间的热量传递，分别计算各个节点的总热流量；计算任意一节点的总热流量时，断裂单元和断裂边的换热系数需要乘以一个折减系数；根据相邻实体单元间的温度差和断裂单元或断裂边的乘以折减系数的换热系数，计算得到实体单元间经过所述任意一节点传递的热流量；根据实体单元的所有节点的温度，求得实体单元内的温度梯度，然后求得实体单元内流入所述任意一节点的热流量；将上述经过所述任意一节点传递的热流量和流入所述任意一节点的热流量相加得到所述任意一节点的总热流量。以三角形单元离散为例来予以说明，如图5所示，用三角形单元的3个节点的温度来代表该三角形单元的温度；于是，整个连续体的温度由所有三角形单元的节点的温度来表征。基于图5所示的拓扑连接图来计算整个连续体的温度场的过程如下。

如图5所示，将三角形单元的节点分别标记为节点1、节点2和节点3等，以图5中的节点1为例，节点1所在的三角形单元有Δ123与三角单元Δ456、Δ789通过节理单元相连接。由于节点1的温度与该三角形单元的节点2、3的温度可能不同，因此在该三角形单元内可能会发生热传导。另外，节点1的温度还可能与节点7、4的温度不同，因此节点1和节点7、4还可能通过节理单元发生热量传递。以与节点1直接相连的三角形单元Δ123为例，设这三个节点的温度为T_i，节点坐标为(x_i,y_i)，i为节点号，i＝1,2,3。假设三角形单元内温度场分布服从线性分布，那么同一三角形单元内任意一点的温度梯度为常量，可表示为：

根据高斯散度定理，公式(1)可写为：

其中，A是三角形单元的面积，n_i是外法向向量，

是m边的平均温度，

是m边的两顶点的坐标分量之差，∈_ij是二维置换张量，

将式(2)代入热传导傅里叶定律，

即可得到沿x,y方向的热流速q_x和q_y；

于是，单位时间内，流入节点1的热流量通过公式(3)计算：

其中

为节点n所对边的单位外法向向量，L⁽ⁿ⁾为节点n所对边的边长，这样便求得了三角形单元Δ123中流入节点1的热流量Q_Δ123→1。

另外，由于三角形单元Δ123与Δ456、Δ789通过节理单元相连接，Δ123与Δ456、Δ789会通过节理单元发生热量交换。以三角形单元Δ123与Δ456为例，绘制如图6所示的示意图，设节点1、2、4、6的温度为T₁,T₂,T₄,T₆，节理单元1264的连接的两个三角形单元Δ123与Δ456，该节理单元的换热系数为h_j，于是单位时间内，三角形单元Δ123与Δ456通过节理单元1264交换的热量为：

其中，L为节理单元的边长。

于是，从三角形单元Δ456通过节理单元流入Δ123的热流量分配到节点1的值为：

与(5)式类似，我们可求得从三角形单元Δ789通过节理单元流入Δ123的热流量分配到节点1的值为：

于是，单位时间内，流入节点1的总的热流量为

Q_total＝Q_Δ123→1+Q_Δ456→1+Q_Δ789→1 (7)

其余节点的总的热流量可根据上述节点1的总的热流量计算方法类似求得。

S107：根据有限差分积分策略，利用计算得到的所述各个节点的总热流量分别不断更新对应节点的温度，最终完成连续-非连续介质区域的热传导；所述任意一节点的温度可按照以下公式不断进行更新：

其中，T^t+Δt是所述任意一节点在下一个时间步的温度，T^t是所述任意一节点在当前时间步t的温度，Δt是时间步长，C_p是比热容，M是所述任意一节点的质量。

一种连续-非连续介质热传导的二维数值模拟系统，其特征在于：包括网格划分模块、判断模块、插入标记模块、非插入标记模块、计算模块和更新模块；

网格划分模块，用于对连续-非连续介质区域进行网格划分；网格划分所形成的网格为实体单元，网格的各个顶点作为节点；

判断模块，用于判断是否在网格间的公共边上插入了节理单元或者界面单元？

插入标记模块，用于将插入在包括节理和裂隙非连续面处公共边的节理单元或者界面单元标记为断裂单元，将插入在节理和裂隙连续面处公共边的节理单元或者界面单元标记为非断裂单元；

非插入标记模块，用于将位于包括节理和裂隙非连续面处的公共边标记为断裂边，将连续面处的公共边标记为非断裂边；

计算模块，用于分别计算各个节点的总热流量；计算任意一节点的总热流量时，所述断裂单元和断裂边的换热系数需要乘以一个折减系数；根据相邻实体单元间的温度差和所述断裂单元或断裂边的乘以折减系数的换热系数，计算得到实体单元间经过所述任意一节点传递的热流量；根据实体单元的所有节点的温度，求得实体单元内的温度梯度，然后求得实体单元内流入所述任意一节点的热流量；将上述经过所述任意一节点传递的热流量和流入所述任意一节点的热流量相加得到所述任意一节点的总热流量；

更新模块，用于计算得到的所述各个节点的总热流量分别不断更新对应节点的温度，最终完成连续-非连续介质区域的热传导；所述任意一节点的温度可按照以下公式不断进行更新：

其中，T^t+Δt是所述任意一节点在下一个时间步的温度，T^t是所述任意一节点在当前时间步t的温度，Δt是时间步长，C_p是比热容，M是所述任意一节点的质量。

本发明的有益效果是：本发明将连续-非连续介质进行网格区域化为实体单元，使计算过程更为容易理解、也更简单；通过乘以折减系数，考虑到裂缝对传热的阻碍作用和在裂缝的动态扩展过程中裂缝对热传导的影响，结合FDEM、DEM、GDEM、DDA、NMM、FEM结合构建热力耦合模型，用以模拟热破裂，提高了计算结果的准确性，本发明提供的技术方案操作简单，适用性和实用性强。

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

Claims (5)

1.一种连续-非连续介质热传导的二维数值模拟方法，其特征在于：包括以下步骤：

S101：对连续-非连续介质区域进行网格划分；网格划分所形成的网格为实体单元，网格的各个顶点作为节点；

S102：判断是否在网格间的公共边上插入了节理单元或者界面单元？若是，则到步骤S103，若否，则到步骤S104；公共边分为节理和裂隙非连续面处公共边和连续面处公共边；

S103：将插入在包括节理和裂隙非连续面处公共边的节理单元或者界面单元标记为断裂单元，将插入在连续面处公共边的节理单元或者界面单元标记为非断裂单元，然后到步骤S105；

S104：将位于包括节理和裂隙非连续面处的公共边标记为断裂边，将连续面处的公共边标记为非断裂边，然后到步骤S105；

S105：通过所述断裂单元或者非断裂单元、断裂边或者非断裂边，实体单元内发生热传导，相邻实体单元间发生热量传递；

S106：根据实体单元内的热传导和相邻实体单元间的热量传递，分别计算各个节点的总热流量；

S107：根据有限差分积分策略，利用计算得到的所述各个节点的总热流量分别不断更新对应节点的温度，最终完成连续-非连续介质区域的热传导。

2.如权利要求1所述的一种连续-非连续介质热传导的二维数值模拟方法，其特征在于：在步骤S106中，计算任意一节点的总热流量时，断裂单元和断裂边的换热系数需要乘以一个折减系数；根据相邻实体单元间的温度差和断裂单元或断裂边的乘以折减系数的换热系数，计算得到实体单元间经过所述任意一节点传递的热流量；根据实体单元的所有节点的温度，求得实体单元内的温度梯度，然后求得实体单元内流入所述任意一节点的热流量；将上述经过所述任意一节点传递的热流量和流入所述任意一节点的热流量相加得到所述任意一节点的总热流量。

3.如权利要求1所述的一种连续-非连续介质热传导的二维数值模拟方法，其特征在于：在步骤S107中，所述任意一节点的温度可按照以下公式不断进行更新：

其中，T^t+Δt是所述任意一节点在下一个时间步的温度，T^t是所述任意一节点在当前时间步t的温度，Δt是时间步长，C_p是比热容，M是所述任意一节点的质量，Q_total是流入节点的热流量。

4.一种连续-非连续介质热传导的二维数值模拟系统，其特征在于：包括网格划分模块、判断模块、插入标记模块、非插入标记模块、计算模块和更新模块；

网格划分模块，用于对连续-非连续介质区域进行网格划分；网格划分所形成的网格为实体单元，网格的各个顶点作为节点；

判断模块，用于判断是否在网格间的公共边上插入了节理单元或者界面单元？

插入标记模块，用于将插入在包括节理和裂隙非连续面处公共边的节理单元或者界面单元标记为断裂单元，将插入在连续面处公共边的节理单元或者界面单元标记为非断裂单元；

非插入标记模块，用于将位于包括节理和裂隙非连续面处的公共边标记为断裂边，将连续面处的公共边标记为非断裂边；

计算模块，用于分别计算各个节点的总热流量；计算任意一节点的总热流量时，所述断裂单元和断裂边的换热系数需要乘以一个折减系数；根据相邻实体单元间的温度差和所述断裂单元或断裂边的乘以折减系数的换热系数，计算得到实体单元间经过所述任意一节点传递的热流量；根据实体单元的所有节点的温度，求得实体单元内的温度梯度，然后求得实体单元内流入所述任意一节点的热流量；将上述经过所述任意一节点传递的热流量和流入所述任意一节点的热流量相加得到所述任意一节点的总热流量；

更新模块，用于计算得到的所述各个节点的总热流量分别不断更新对应节点的温度，最终完成连续-非连续介质区域的热传导。

5.如权利要求4所述的一种连续-非连续介质热传导的二维数值模拟系统，其特征在于：在更新模块中，所述任意一节点的温度可按照以下公式不断进行更新：

其中，T^t+Δt是所述任意一节点在下一个时间步的温度，T^t是所述任意一节点在当前时间步t的温度，Δt是时间步长，C_p是比热容，M是所述任意一节点的质量，Q_total是流入节点的热流量。

Images