计算热物理

目录题目 (2)1.方程和定解条件的无量纲化 (4)2.控制方程的离散 (4)3.稳定性分析 (5)4.温度场的计算求解 (6)5.温度场的等值线图 (9)6.分析比较 (12)7.个人心得体会 (14)附录 (16)题目：定义在正方形区域（如图所示）H*H=3m*3m 的二维非稳态导热问题的控制方程及其定解条件为其中ρ，c，λ分别为导热体的密度、比热和导热系数，且ρ=7820kg/m3，c=460J/(kg•K)，λ=15W/(m•K)。

将定义域各方向均匀分成10个区块，取点中心网格分割方式，各边界共11个节点; 或者取块中心网格分割，各边界共10个节点，（两种网格，至少做一种），按以下要求，数值求解该问题。

（1）取（X，Y）=（x，y）H ，θ=T−T0T1−T0，τ=tρcH2/λ，将方程和定解条件无量纲化，并在无量纲化条件下求解。

（2）采用以下方法对控制方程离散：（a）用控制容积法离散成显式、全隐式格式；（b）用ADI格式离散。

（3）用Fourier分析方法对以上3种格式的稳定性进行分析。

（4）嵌入初、边条件，分别对以上3种格式构成的定解问题编制计算程序，作数值求解，直到达到稳态温度分布——能直接求解显式格式和ADI格式构成的代数方程）直接求解，不能直接求解的（隐式格式构成的代数方程）采用以下方法迭代求解：（a）Jacobi简单点迭代；（b）Guess-Seidel点迭代；（c）带松弛的Guess-Seidel点迭代；（d）带松弛的Guess-Seidel线迭代；（e）基于Guess-Seidel的交替方向线迭代。

并对收敛结果做出比较。

（5）选择全隐式格式在5个不同时刻（包括1个到达稳态的时刻）的计算结果，将无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。

从5个不同时刻的等值线图，分析时间发展中T的变化特性。

（6）根据数值计算情况，（精度高低，CPU时间消耗，编制程序难易等）对各种不同的离散方法及代数解法进行分析比较，并做出结论。

（7）综合上述内容和作业情况，写出作业报告，并对此类作业提出建议和要求。

解：（1）方程和定解条件的无量纲化：将（X ，Y ）=（x ，y ）H，θ=T−T 0T 1−T 0，τ=tρcH 2/λ代入原方程组，无量纲化，得到：（2） 控制方程的离散： A ． 控制容积法假定非稳态项中温度随空间为阶梯分布，扩散项中温度随时间线性分布、随空间阶梯分布 显式格式:(θP −θP 0)∆X∆Y=(θE 0−θP(δX )e−θP 0−θW(δX )w)∆Y∆τ+(θN 0−θP 0(δY )n−θP 0−θS(δY )s)∆X∆τ写成：a P θP =a E θE 0+a W θW 0+a N θN 0+a S θS 0+b 的形式，则：a E =∆Y(δX )e ，a W =∆Y(δX )w，a N =∆X(δY )n，a S =∆X(δY )sa P =∆X∆Y∆τ，a P 0=∆X∆Y ∆τ−a E −a W −a N −a S ，b =a P 0θP 0全隐式格式：(θP −θP 0)∆X∆Y =(θE−θP (δX )e−θP −θW (δX )w)∆Y∆τ+(θN−θP(δY )n−θP −θS (δY )s)∆X∆τ写成：a P θP =a E θE +a W θW +a N θN +a S θS +b 的形式，则：a E=∆Y(δX)e ，a W=∆Y(δX)w，a N=∆X(δY)n，a S=∆X(δY)sa P=∆X∆Y∆τ+a E+a W+a N+a S，a P0=∆X∆Y∆τ，b=a P0θP0B．ADI格式：θi,j∗−θi,j n ∆τ/2=θi+1,j∗−2θi,j∗+θi−1,j∗∆X2+θi,j+1n−2θi,jn+θi,j−1n∆Y2θi,j n+1−θi,j∗∆τ/2=θi+1,j∗−2θi,j∗+θi−1,j∗∆X2+θi,j+1n+1−2θi,jn+1+θi,j−1n+1∆Y2（3）稳定性分析:A．显示格式：a PθP=a EθE0+a WθW0+a NθN0+a SθS0+ba P g=a E e ik X∆X+a W e−ik X∆X+a N e ik Y∆Y+a S e−ik Y∆Y+a P0(*)本题中定义域各方向均匀分成10个区块，则(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1a E=a W=a N=a S=1，a P=∆X∆Y∆τ则(*)式可化为：g=2cos(k X∆X)+2cos(k Y∆Y)+(∆X∆Y/∆τ−4)∆X∆Y/∆τ假设k X=k Y，则g=1−8∆τ∆X2sin2(k∆X2)有|g|≤1恒成立，则∆τ∆X2sin2(k∆X2)≤14∆τ∆X2≤14，∆τ≤14∆X2。

即∆τ≤14∆X2时，该显式格式收敛。

本题中，取∆X=0.1，故∆τ≤0.0025。

B．全隐式格式：a PθP=a EθE+a WθW+a NθN+a SθS+ba P g=g(a E e ik X∆X+a W e−ik X∆X+a N e ik Y∆Y+a S e−ik Y∆Y)+a P0(**)本题中定义域各方向均匀分成10个区块，则(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1a E=a W=a N=a S=1，a P0=∆X∆Y∆τ(**)式化为：g=∆X∆Y/∆τ4−2cos(k X∆X)−2cos(k Y∆Y)+∆X∆Y/∆τ显然|g|≤1恒成立，故全隐式格式恒稳定。

C．ADI格式：θi,j∗−θi,j n ∆τ/2=θi+1,j∗−2θi,j∗+θi−1,j∗∆X2+θi,j+1n−2θi,jn+θi,j−1n∆Y2θi,j n+1−θi,j∗∆τ/2=θi+1,j∗−2θi,j∗+θi−1,j∗∆X2+θi,j+1n+1−2θi,jn+1+θi,j−1n+1∆Y2空间导数离散式用微分算子符号写出，上面两个等式可化为：(1−αX2δX2)θi,j∗=(1+αY2δY2)θi,j n，αX=∆τ/∆X2(1−αY2δY2)θi,j n+1=(1+αX2δX2)θi,j∗，αY=∆τ/∆Y2消去θi,j∗，得：(1−αX2δX2)(1−αY2δY2)θi,j n+1=(1+αX2δX2)(1+αY2δY2)θi,j n故g=(1−2αX sin2k X∆X2)(1−2αY sin2k Y∆Y2)(1+2αX sin2k X∆X2)(1+2αY sin2k Y∆Y2)显然|g|≤1恒成立，故ADI格式恒稳定。

（4）温度场的计算求解将定义域各方向均匀分成10 个区块,采用点中心网格分割方式，各边界共11 个节点。

边界条件的离散：（1）Y=0：θi,1=1(i=1—11)；（2）Y=1：θi,11=0(i=1—11)；（3）X=0：j=2—10显式：a 1,j θ1,j =a 2,j θ2,j +a 1,j+1θ1,j+1+a 1,j−1θ1,j−1+b 其中： a 1,j =∆X∆Y 2∆τ，a 2,j =∆Y(δX )e，a 1,j+1=∆X 2(δY )n，a 1,j−1=∆X 2(δY )sb =θ1,j 0(a 1,j −a 2,j −a 1,j+1−a 1,j−1)+q w ∆Y隐式：a 1,j θ1,j =a 2,j θ2,j +a 1,j+1θ1,j+1+a 1,j−1θ1,j−1+b 其中： a 2,j =∆Y (δX )e ，a 1,j+1=∆X 2(δY )n，a 1,j−1=∆X 2(δY )sa 1,j =∆X∆Y2∆τ+ a 2,j +a 1,j+1+a 1,j−1，b =∆X∆Y 2∆τθ1,j 0+q w ∆Y（4）X=1：j=2—10显式：a 11,j θ11,j =a 10,j θ10,j +a 11,j+1θ11,j+1+a 11,j−1θ11,j−1+b 其中：a 11,j =∆X∆Y 2∆τ，a 10,j =∆Y(δX )w，a 11,j+1=∆X 2(δY )n，a 11,j−1=∆X 2(δY )sb =θ11,j 0(a 11,j −a 10,j −a 11,j+1−a 11,j−1)+q e ∆Y隐式：a 11,j θ11,j =a 10,j θ10,j +a 11,j+1θ11,j+1+a 11,j−1θ11,j−1+b 其中： a 10,j =∆Y(δX )w，a 11,j+1=∆X 2(δY )n，a 11,j−1=∆X 2(δY )sa 11,j =∆X∆Y 2∆τ+ a 10,j +a 11,j+1+a 11,j−1，b =∆X∆Y 2∆τθ11,j 0+q e ∆Y初始条件：由“初始T(x,y,0)为x=0→x=H 的温度成线性分布”我们可以推知：θi,j =1−0.1∗(i −1) i,j=1—11据此编程，编程软件：Visual C++，程序代码见附录。

其中，左右边界的热流量qe 和qw 、时间步长、稳态精度、迭代精度、时间上限以及最大迭代次数均可在define 项目中直接调整。

几种迭代算法的比较：1.Jacobi简单点迭代：利用上一时刻的节点温度值计算下一时刻该节点的温度，程序较为简单。

收敛速度慢，迭代最大次数为9次，出现在时间刚开始的阶段。

下面是计算结果：2.Gauss-Seidel点迭代：在计算节点值时，采用邻近节点的最新可用值，而不局限于前一迭代完成后所得到的值。

此种方法由于及时引入新值，迭代收敛速度快于简单迭代。

在实际程序设计中，先计算边界，使边界的最新可用值被尽可能快的引入到内节点的计算中，对加快收敛速度起到很好的作用。

在以后的程序中，也都采用了这种方法。

该算法的最大迭代次数为7次，出现在时间值很小的起步阶段。

下面是计算结果：3.松弛迭代本次上机取松弛因子w在从0.1开始按0.1增量递增到1.9。

松弛迭代的收敛速度与松弛因子有关。

当松弛因子在1.2-1.4时，收敛速度明显加快，而当松弛因子取0.1或1.9附近的值时，有可能发散。

4.Gauss-Seidel的交替方向线迭代此法扫描方向可以变化，且有多种组合。

本程序先处理边界条件，然后采用从左至右的列扫描，后再采用从上至下的行扫描，全场两次扫描组成一轮迭代，这种扫描是对一条线上隐式格式迭代求解。

相对于Gauss-Seidel迭代，收敛速度有所提高，最大迭代次数为4次。

下面是计算结果：（5）温度场的等值线图这里，结合Gauss-Seidel的交替方向线迭代的计算结果，无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。

我们可以看到，达到稳态所需的无量纲时间为0.636，折算成实际时间为38小时。