基本原理

对于多维非稳态导热问题，其时间项的离散格式有两种：显式格式和隐式格式。在显式格式中，下一时间步的温度场由上一时间步的已知量求出，不需要求解代数方程组，但其时间步的选取受稳定性条件的限制，故通常只能取很小的时间步。隐式格式不受稳定性限制，但其需要对五对角矩阵（二维）或七对角矩阵（三维）进行求解，故每一时间步的计算量较大。

相比之下，交替方向隐式方法^[1]（alternating direction implicit, ADI）是在显式与隐式之间折中的一种方法。对于二维问题，先将xxx方向的温度值按隐式处理，yyy方向的温度值按显式处理，这样就将二维问题转化为一系列并列的一维问题，可以直接使用三对角阵算法（TDMA）进行计算；然后再倒过来，将yyy方向隐式，xxx方向显式处理。这种方法被称为Peaceman-Rachford ADI格式。可以证明，该方法对于二维问题是绝对稳定的，故没有稳定性限制，可以取较大的时间步长，同时每一时间步的计算量较小，可大幅加快计算速度。

下面介绍Peaceman-Rachford ADI格式下二维非稳态导热方程的求解过程。

基于有限体积法的全隐离散方程可表示为：
aP0(TP−TP0)=aE(TE−TP)−aW(TP−TW)+aN(TN−TP)−aS(TP−TS)+SΔxΔya_P^0(T_P-T_P^0)=a_E(T_E-T_P)-a_W(T_P-T_W)+a_N(T_N-T_P)-a_S(T_P-T_S)+S\Delta x\Delta y aP0(TP−TP0)=aE(TE−TP)−aW(TP−TW)+aN(TN−TP)−aS(TP−TS)+SΔxΔy

(1)沿xxx方向求解：

先取Δt/2\Delta t/2Δt/2的时间步长，对x方向做隐式处理，求解温度场的中间量UUU，则全隐离散方程可改成：
aP0(UP−TP0)=aE(UE−UP)−aW(UP−UW)+aN(TN0−TP0)−aS(TP0−TS0)+SΔxΔya_P^0(U_P-T_P^0)=a_E(U_E-U_P)-a_W(U_P-U_W)+a_N(T_N^0-T_P^0)-a_S(T_P^0-T_S^0)+S\Delta x\Delta y aP0(UP−TP0)=aE(UE−UP)−aW(UP−UW)+aN(TN0−TP0)−aS(TP0−TS0)+SΔxΔy

整理后得到：

(aP0+aE+aW)UP=aEUE+aWUW+(aP0−aN−aS)TP0+aNTN0+aSTS0+SΔxΔy(a_P^0+a_E+a_W)U_P=a_EU_E+a_WU_W+(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y (aP0+aE+aW)UP=aEUE+aWUW+(aP0−aN−aS)TP0+aNTN0+aSTS0+SΔxΔy

在上式中，除了UPU_PUP，UEU_EUE和UWU_WUW，其他量都是已知的。令aP=aP0+aE+aWa_P=a_P^0+a_E+a_WaP=aP0+aE+aW，D=(aP0−aN−aS)TP0+aNTN0+aSTS0+SΔxΔyD=(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta yD=(aP0−aN−aS)TP0+aNTN0+aSTS0+SΔxΔy，则上式可以简化为：

aPUP=aEUE+aWUW+Da_PU_P=a_EU_E+a_WU_W+D aPUP=aEUE+aWUW+D

通过使用沿xxx方向的TDMA可以求得温度场中间量UUU。

(2)沿yyy方向求解：

在求得中间量UUU的基础上，再取Δt/2\Delta t/2Δt/2的时间步长，对y方向做隐式处理，求解下一时刻温度场T1T^1T1，则全隐离散方程可改成：
aP0(TP−TP0)=aE(TE1−TP1)−aW(TP1−TW1)+aN(UN−UP)−aS(UP−US)+SΔxΔya_P^0(T_P-T_P^0)=a_E(T_E^1-T_P^1)-a_W(T_P^1-T_W^1)+a_N(U_N-U_P)-a_S(U_P-U_S)+S\Delta x\Delta y aP0(TP−TP0)=aE(TE1−TP1)−aW(TP1−TW1)+aN(UN−UP)−aS(UP−US)+SΔxΔy

整理后得到：

(aP0+aN+aS)TP1=aNTN1+aSUS1+(aP0−aE−aW)UP+aEUE+aEUW+SΔxΔy(a_P^0+a_N+a_S)T_P^1=a_NT_N^1+a_SU_S^1+(a_P^0-a_E-a_W)U_P+a_EU_E+a_EU_W+S\Delta x\Delta y (aP0+aN+aS)TP1=aNTN1+aSUS1+(aP0−aE−aW)UP+aEUE+aEUW+SΔxΔy

在上式中，除了Tp1T_p^1Tp1，TN1T_N^1TN1和TS1T_S^1TS1，其他量都是已知的。令aP=aP0+aN+aSa_P=a_P^0+a_N+a_SaP=aP0+aN+aS，D=(aP0−aE−aW)UP+aEUE+aEUW+SΔxΔyD=(a_P^0-a_E-a_W)U_P+a_EU_E+a_EU_W+S\Delta x\Delta yD=(aP0−aE−aW)UP+aEUE+aEUW+SΔxΔy，则上式可以简化为：

aPTP1=aNTN1+aNTN1+Da_PT_P^1=a_NT_N^1+a_NT_N^1+D aPTP1=aNTN1+aNTN1+D

再使用沿yyy方向的TDMA求得下一时刻的温度场T1T^1T1。

边界节点的处理方式与上述过程相似，主要区别在于该节点的离散方程中未知项少了一项，例如对于顶部对流边界，aNa_NaN的取值得改为αΔx\alpha\Delta xαΔx，TNT_NTN由未知的待求解量改为已知量TfT_fTf。

求解程序

下面设计一个算例，并通过Fortran程序给出其求解过程：

有一个宽400mm，高200mm的矩形计算域。计算域内材料的密度为2500kg/m³，热容为1372J/(kg·K)。含内热源，热生成率恒定为1050J/(m³·s)。导热系数为温度的多项式函数：

λ=−1.13588×10−15T5+3.25358×10−12T4−3.25305×10−9T3+1.32926×10−6T2−9.27637×10−5T+1.04478\lambda=-1.13588\times 10^{-15}T^5+3.25358\times 10^{-12}T^4-3.25305\times 10^{-9}T^3+1.32926\times 10^{-6}T^2-9.27637\times 10^{-5}T+1.04478 λ=−1.13588×10−15T5+3.25358×10−12T4−3.25305×10−9T3+1.32926×10−6T2−9.27637×10−5T+1.04478

计算域的初始温度为1500℃，四周与环境进行对流换热，其中顶部对流给热系数α1\alpha_1α1为200W/(m²·K)，底部和左右侧的对流给热系数α2\alpha_2α2为80W/(m²·K)，环境温度TfT_fTf为20℃。模拟的物理时间长度为3h。

采用均匀结构化网格，每个网格节点的尺寸为2mm×2mm。由于该算例的模型沿yyy轴是对称的，实际计算域可以只取左半部分，即选取的计算域尺寸为200mm×200mm，网格节点数量为101×101。计算域右侧作为对称边界处理，该边界上TE=TWT_E=T_WTE=TW，aE=aWa_E=a_WaE=aW。

由于Peaceman-Rachford ADI格式的时间步长选取不受稳定性限制，其时间步长的选取比较灵活。由于非稳态导热随着时间的进行，会趋向于稳态，刚开始时温度变化率会很大，随后则会大幅减小，本算例选取时间步长的准则是使选定的时间步长内最大温度变化幅度不超过10℃，这样刚开始时间步长会很小，随后则会变得很大。

计算过程中选取了6个采样点，并设定每隔5分钟输出一次这6个采样点的温度值。

求解程序如下，模拟用到的参数在网格、材料、边界条件、初始条件、求解控制模块中通过常量定义出来，以便于参数的调整：

module meshimplicit nonereal, parameter :: width = 200e-3, height = 200e-3real, parameter :: delta_X = 2e-3, delta_Y = 2e-3
end module meshmodule materialimplicit nonereal, parameter :: rho = 2500, Cp = 1372, S = 1050
end module materialmodule boundary_conditionimplicit nonereal, parameter :: alpha1 = 300, alpha2 = 80real, parameter :: Tf = 20
end module boundary_conditionmodule initial_conditionimplicit nonereal, parameter :: T0 = 1500
end module initial_conditionmodule solution_controlimplicit nonereal, parameter :: time_length = 180 * 60real, parameter :: sampling_interval = 5 * 60real, parameter :: max_deltaT_per_step = 10real, dimension(6) :: sample_points_X = (/30e-3, 100e-3, 200e-3, 30e-3, 100e-3, 200e-3/)real, dimension(6) :: sample_points_Y = (/100e-3, 100e-3, 100e-3, 170e-3, 170e-3, 170e-3/)integer :: spIndex(6,2)
end module solution_controlprogram TDMAuse meshuse materialuse solution_controlimplicit noneinteger :: Nx, Ny!确定X、Y方向的网格节点数量Nx = width / delta_X + 1Ny = height / delta_Y + 1!根据采样点位置获得采样点在网格节点中的索引spIndex(:,1) = floor(sample_points_X / delta_X + 1)spIndex(:,2) = floor(sample_points_Y / delta_Y + 1)!开始计算call calculate(Nx, Ny)end program TDMAsubroutine calculate(Nx, Ny)use meshuse materialuse boundary_conditionuse initial_conditionuse solution_controlimplicit noneinteger, intent(in) :: Nx, Nyinteger :: ireal :: current_time, current_sampling_time, delta_Taureal, dimension(0:Nx+1,0:Ny+1) :: T, U, T1, delta_Treal, dimension(Nx,Ny) :: lambda, aP0, aE, aW, aN, aS, aP, Sv, D, P, Qreal :: kx(Ny), ky(Nx)real :: max_delta_T, max_ratereal :: tStart, tEnd!初始化系数矩阵AaP0 = 0aE = 0aW = 0aN = 0aS = 0!初始化温度矩阵T = TfT(1:Nx,1:Ny) = T0U = TfT1 = Tf!初始化源项矩阵Sv = S * delta_X * delta_YSv(1,:) = Sv(1,:) / 2Sv(:,1) = Sv(:,1) / 2Sv(:,Ny) = Sv(:,Ny) / 2!初始化输出write(*,"(A4,$)") "Time"do i = 1, 5write(*,"(A9,I1,$)") "SP", iend dowrite(*,"(A10)") "SP6"!按时间步进行温度场计算current_time = 0current_sampling_time = sampling_intervaldelta_Tau = 0.0001  !考虑到降温刚开始时的温度变化率很大，初始delta_Tau取一个很小的值call CPU_TIME(tStart)do while (current_time < time_length)!更新当前时间current_time = current_time + delta_Tau!通过多项式函数计算材料导热系数lambda = -1.13588E-15 * T(1:Nx,1:Ny)**5 + 3.25358E-12 * T(1:Nx,1:Ny)**4 - 3.25305E-09 * T(1:Nx,1:Ny)**3 &+ 1.32926E-06 * T(1:Nx,1:Ny)**2 - 9.27637E-05 * T(1:Nx,1:Ny) + 1.04478!计算系数矩阵A!(1)内部节点的计算aP0 = rho * Cp * delta_X * delta_Y / (delta_Tau / 2)aE(1:Nx-1,:) = 2 / (1 / lambda(1:Nx-1,:) + 1 / lambda(2:Nx,:)) / delta_X * delta_YaW(2:Nx,:) = 2 / (1 / lambda(2:Nx,:) + 1 / lambda(1:Nx-1,:)) / delta_X * delta_YaN(:,1:Ny-1) = 2 / (1 / lambda(:,1:Ny-1) + 1 / lambda(:,2:Ny)) / delta_Y * delta_XaS(:,2:Ny) = 2 / (1 / lambda(:,2:Ny) + 1 / lambda(:,1:Ny-1)) / delta_Y * delta_X!(2)左侧、上下侧对流边界节点的处理aW(1,:) = alpha2 * delta_YaN(:,Ny) = alpha1 * delta_XaS(:,1) = alpha2 * delta_X!(3)右侧对称边界的处理aW(Nx,:) = aW(Nx,:) * 2!(4)边界节点的减半处理!(4-1)上侧aP0、aE、aW减半aP0(:,Ny) = aP0(:,Ny) / 2aE(:,Ny) = aE(:,Ny) / 2aW(:,Ny) = aW(:,Ny) / 2!(4-2)下侧aP0、aE、aW减半aP0(:,1) = aP0(:,1) / 2aE(:,1) = aE(:,1) / 2aW(:,1) = aW(:,1) / 2!(4-3)左侧aP0、aN、aS减半aP0(1,:) = aP0(1,:) / 2aN(1,:) = aN(1,:) / 2aS(1,:) = aS(1,:) / 2!沿x方向使用TDMA隐式求解aP = aP0 + aE + aW!(1)定义TDMA中的常数矩阵D = (aP0 - aN - aS) * T(1:Nx,1:Ny) + aN * T(1:Nx,2:Ny+1) + aS * T(1:Nx,0:Ny-1) + Sv!(2)左侧边界的特殊处理D(1,:) = D(1,:) + aW(1,:) * Tf!!!!!!P(1,:) = aE(1,:) / aP(1,:)Q(1,:) = D(1,:) / aP(1,:)do i = 2, Nxkx = aP(i,:) - aW(i,:) * P(i-1,:)P(i,:) = aE(i,:) / kxQ(i,:) = (D(i,:) + aW(i,:) * Q(i-1,:)) / kxend doU(Nx,1:Ny) = Q(Nx,:)do i = Nx-1, 1, -1U(i,1:Ny) = P(i,:) * U(i+1,1:Ny) + Q(i,:)end do!沿y方向使用TDMA隐式求解aP = aP0 + aN + aS!(1)定义TDMA中的常数矩阵D = (aP0 - aE - aW) * U(1:Nx,1:Ny) + aE * U(2:Nx+1,1:Ny) + aW * U(0:Nx-1,1:Ny) + Sv!(2)上下侧边界的特殊处理D(:,1) = D(:,1) + aS(:,1) * TfD(:,Ny) = D(:,Ny) + aN(:,Ny) * Tf!!!!!!P(:,1) = aN(:,1) / aP(:,1)Q(:,1) = D(:,1) / aP(:,1)do i = 2, Nyky = aP(:,i) - aS(:,i) * P(:,i-1)P(:,i) = aN(:,i) / kyQ(:,i) = (D(:,i) + aS(:,i) * Q(:,i-1)) / kyend doT1(1:Nx,Ny) = Q(:,Ny)do i = Ny-1, 1, -1T1(1:Nx,i) = P(:,i) * T1(1:Nx,i+1) + Q(:,i)end do!当达到采样时间时，输出采样点温度if (current_time >= current_sampling_time) then!输出6个采样点的温度值write(*,"(I4,$)") floor(current_time / 60)do i = 1, 5write(*,"(f10.2,$)") T1(spIndex(i,1), spIndex(i,2))end dowrite(*,"(f10.2)") T1(spIndex(6,1), spIndex(6,2))!!更新下一次需要采样的时间current_sampling_time = current_sampling_time + sampling_intervalend if!计算最大冷却速率delta_T = abs(T1 - T)max_delta_T = maxval(delta_T)max_rate = max_delta_T / delta_Tau!根据对每个时间步内的最大温度变化率限制，以及下一次采样时间确定下一个时间步长delta_Tau = min(max_deltaT_per_step / max_rate, current_sampling_time - current_time)!更新温度矩阵TT = T1end docall CPU_TIME(tEnd)write(*,"('Calculate complete, time consumption: ', f8.4, ' sec')") (tEnd - tStart)end subroutine calculate

参考文献

[1] 陶文铨. 数值传热学[M]. 第二版. 西安: 西安交通大学出版社, 2001: 99-104.

二维导热问题的ADI-TDMA算法相关推荐

二维光子晶体禁带的遗传优化算法实现
二维光子晶体禁带的遗传优化算法MATLAB源代码光子晶体中因周期性结构而存在的频率禁带称为光子禁带,光子禁带的存在是光子晶体具有广泛应用前景的重要原因. 禁带越大,可控光的频带也越宽,因此如何设计合 ...
看懂二维码识别OCR：从算法到API 接入代码
引言二维码识别OCR(Optical Character Recognition)是结合了图像处理和OCR技术,以识别和提取二维码中的信息的技术,二维码识别OCR 可以实现对图像中的二维码进行自动检 ...
基于Excel的QR二维码生成工具——原理及算法详解(之一）
老虎二维码(下载链接在这里)是一个基于Excel的二维码生成工具,完全使用Excel表单公式结合VBA实现,没有调用任何外部库,实现了支持中文英文混合字符以及常用微信二维码编码的自动生成,在工作表单元 ...
数组（一维数组、多维数组/二维数组）和简单排序算法
提示:数组是线性数据结构中最为基础,最为典型的一种顺序型结构. 它用一组连续的内存空间 ,来存储一组具有相同类型的数据. 与变量相比,变量是一种单一的数据存储方式,而数组是用于存储一连串的一组数据. ...
矩阵（二维数组）的性质在算法求解中的应用
本文所说的矩阵(matrix),其实在编程实现时,往往以二维数组的形式出现. 1. 对称矩阵(二维数组) 在求解旅行商问题时,题干中要求,城市之间彼此互通(两城市之间的道路只有一条). double ...
二维点云拉普拉斯深度平滑算法-matlab
克服普通平滑会破坏括扑关系的问题,引入权值,更加可靠输入:data:数据集 kkk:邻域搜索范围,视情况定 M:迭代次数,越高,平滑越夸张 function [data]=Laplacian_dee ...
【数据挖掘】K-Means 二维数据聚类分析 ( K-Means 迭代总结 | K-Means 初始中心点选择方案 | K-Means 算法优缺点 | K-Means 算法变种 )
文章目录 K-Means 二维数据聚类分析数据样本及聚类要求二维数据曼哈顿距离计算 K-Means 算法步骤第一次迭代 : 步骤 ( 1 ) 中心点初始化第一次迭代 : 步骤 ( 2 ) ...
c语言求解热传导方程,二维稳态导热问题的数值解法.docx
核科学与技术学院 <传热学> 二维稳态导热问题的数值解法作业姓名:罗晓学号: 2014151214 班级:任课教师:李磊,张智刚哈尔滨工程大学核科学与技术学院 2016 年 11 ...
种子点生长算法（上）——二维种子点生长
本文只用作笔记,方便日后查阅. 种子点生长算法上--二维种子点生长下文提到的种子点生长算法,包括泛洪法,扫描线法,区段法三种.文本先从最简单的泛洪法入手介绍种子点生长算法的相关概念.之后进一步讨论了 ...
基于Rplidar二维雷达使用Hector_SLAM算法在ROS中建图
文章目录前言一.ROS分布式通信(配置多机通信) 1.简介 2.步骤 2.1 准备 2.2 修改配置文件 2.3配置主机IP 2.4配置从机IP 二.RPlidar的使用教程 1.创建环境 2.下 ...

二维导热问题的ADI-TDMA算法

二维导热问题的ADI-TDMA算法

基本原理

求解程序

参考文献

二维导热问题的ADI-TDMA算法相关推荐

最新文章

热门文章