计算热物理

目录题目...........................................................................................21.方程和定解条件的无量纲化................................................42.控制方程的离散....................................................................43.稳定性分析............................................................................54.温度场的计算求解................................................................65.温度场的等值线图................................................................96.分析比较..............................................................................127.个人心得体会......................................................................14附录.........................................................................................16计算热物理大作业1题目：定义在正方形区域（如图所示）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（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）综合上述内容和作业情况，写出作业报告，并对此类作业提出计算热物理大作业3建议和要求。解：（1）方程和定解条件的无量纲化：将（X，Y）=（x，y）H，θ=T−T0T1−T0，τ=tρcH2/λ代入原方程组，无量纲化，得到：（2）控制方程的离散：A．控制容积法假定非稳态项中温度随空间为阶梯分布，扩散项中温度随时间线性分布、随空间阶梯分布显式格式:(θP−θP0)∆X∆Y=(θE0−θP0(δX)e−θP0−θW0(δX)w)∆Y∆τ+(θN0−θP0(δY)n−θP0−θS0(δY)s)∆X∆τ写成：aPθP=aEθE0+aWθW0+aNθN0+aSθS0+b的形式，则：aE=∆Y(δX)e，aW=∆Y(δX)w，aN=∆X(δY)n，aS=∆X(δY)saP=∆X∆Y∆τ，aP0=∆X∆Y∆τ−aE−aW−aN−aS，b=aP0θP0全隐式格式：(θP−θP0)∆X∆Y=(θE−θP(δX)e−θP−θW(δX)w)∆Y∆τ+(θN−θP(δY)n−θP−θS(δY)s)∆X∆τ写成：aPθP=aEθE+aWθW+aNθN+aSθS+b的形式，计算热物理大作业4则：aE=∆Y(δX)e，aW=∆Y(δX)w，aN=∆X(δY)n，aS=∆X(δY)saP=∆X∆Y∆τ+aE+aW+aN+aS，aP0=∆X∆Y∆τ，b=aP0θP0B．ADI格式：θi,j∗−θi,jn∆τ/2=θi+1,j∗−2θi,j∗+θi−1,j∗∆X2+θi,j+1n−2θi,jn+θi,j−1n∆Y2θi,jn+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．显示格式：aPθP=aEθE0+aWθW0+aNθN0+aSθS0+baPg=aEeikX∆X+aWe−ikX∆X+aNeikY∆Y+aSe−ikY∆Y+aP0(*)本题中定义域各方向均匀分成10个区块，则(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1aE=aW=aN=aS=1，aP=∆X∆Y∆τ则(*)式可化为：g=2cos(kX∆X)+2cos(kY∆Y)+(∆X∆Y/∆τ−4)∆X∆Y/∆τ假设kX=kY，则g=1−8∆τ∆X2sin2(k∆X2)有|g|≤1恒成立，则∆τ∆X2sin2(k∆X2)≤14∆τ∆X2≤14，∆τ≤14∆X2。即∆τ≤14∆X2时，该显式格式收敛。本题中，取∆X=0.1，故∆τ≤0.0025。B．全隐式格式：aPθP=aEθE+aWθW+aNθN+aSθS+baPg=g(aEeikX∆X+aWe−ikX∆X+aNeikY∆Y+aSe−ikY∆Y)+aP0(**)计算热物理大作业5本题中定义域各方向均匀分成10个区块，则(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1aE=aW=aN=aS=1，aP0=∆X∆Y∆τ(**)式化为：g=∆X∆Y/∆τ4−2cos(kX∆X)−2cos(kY∆Y)+∆X∆Y/∆τ显然|g|≤1恒成立，故全隐式格式恒稳定。C．ADI格式：θi,j∗−θi,jn∆τ/2=θi+1,j∗−2θi,j∗+θi−1,j∗∆X2+θi,j+1n−2θi,jn+θi,j−1n∆Y2θi,jn+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,jn，αX=∆τ/∆X2(1−αY2δY2)θi,jn+1=(1+αX2δX2)θi,j∗，αY=∆τ/∆Y2消去θi,j∗，得：(1−αX2δX2)(1−αY2δY2)θi,jn+1=(1+αX2δX2)(1+αY2δY2)θi,jn故g=(1−2αXsin2kX∆X2)(1−2αYsin2kY∆Y2)(1+2αXsin2kX∆X2)(1+2αYsin2kY∆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)；计算热物理大作业6（3）X=0：j=2—10显式：a1,jθ1,j=a2,jθ2,j+a1,j+1θ1,j+1+a1,j−1θ1,j−1+b其中：a1,j=∆X∆Y2∆τ，a2,j=∆Y(δX)e，a1,j+1=∆X2(δY)n，a1,j−1=∆X2(δY)sb=θ1,j0(a1,j−a2,j−a1,j+1−a1,j−1)+qw∆Y隐式：a1,jθ1,j=a2,jθ2,j+a1,j+1θ1,j+1+a1,j−1θ1,j−1+b其中：a2,j=∆Y(δX)e，a1,j+1=∆X2(δY)n，a1,j−1=∆X2(δY)sa1,j=∆X∆Y2∆τ+a2,j+a1,j+1+a1,j−1，b=∆X∆Y2∆τθ1,j0+qw∆Y（4）X=1：j=2—10显式：a11,jθ11,j=a10,jθ10,j+a11,j+1θ11,j+1+a11,j−1θ11,j−1+b其中：a11,j=∆X∆Y2∆τ，a10,j=∆Y(δX)w，a11,j+1=∆X2(δY)n，a11,j−1=∆X2(δY)sb=θ11,j0(a11,j−a10,j−a11,j+1−a11,j−1)+qe∆Y隐式：a11,jθ11,j=a10,jθ10,j+a11,j+1θ11,j+1+a11,j−1θ11,j−1+b其中：a10,j=∆Y(δX)w，a11,j+1=∆X2(δY)n，a11,j−1=∆X2(δY)sa11,j=∆X∆Y2∆τ+a10,j+a11,j+1+a11,j−1，b=∆X∆Y2∆τθ11,j0+qe∆Y初始条件：由“初始T(x,y,0)为x=0→x=H的温度成线性分布”我们可以推知：θi,j=1−0.1∗(i−1)i,j=1—11据此编程，编程软件：VisualC++，程序代码见附录。其中，左右边界的热流量qe和qw、时间步长、稳态精度、迭代精度、时间上限以及最大迭代次数均可在define项目中直接调整。几种迭代算法的比较：计算热物理大作业71.Jacobi简单点迭代：利用上一时刻的节点温度值计算下一时刻该节点的温度，程序较为简单。收敛速度慢，迭代最大次数为9次，出现在时间刚开始的阶段。下面是计算结果：2.Gauss-Seidel点迭代：在计算节点值时，采用邻近节点的最新可用值，而不局限于前一迭代完成后所得到的值。此种方法由于及时引入新值，迭代收敛速度快于简单迭代。在实际程序设计中，先计算边界，使边界的最新可用值被尽可能快的引入到内节点的计算中，对加快收敛速度起到很好的作用。在以后的程序中，也都采用了这种方法。该算法的最大迭代次数为7次，出现在时间值很小的起步阶段。下面是计算结果：计算热物理大作业83.松弛迭代本次上机取松弛因子w在从0.1开始按0.1增量递增到1.9。松弛迭代的收敛速度与松弛因子有关。当松弛因子在1.2-1.4时，收敛速度明显加快，而当松弛因子取0.1或1.9附近的值时，有可能发散。4.Gauss-Seidel的交替方向线迭代此法扫描方向可以变化，且有多种组合。本程序先处理边界条件，然后采用从左至右的列扫描，后再采用从上至下的行扫描，全场两次扫描组成一轮迭代，这种扫描是对一条线上隐式格式迭代求解。相对于Gauss-Seidel迭代，收敛速度有所提高，最大迭代次数为4次。下面是计算结果：计算热物理大作业9（5）温度场的等值线图这里，结合Gauss-Seidel的交替方向线迭代的计算结果，无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。我们可以看到，达到稳态所需的无量纲时间为0.636，折算成实际时间为38小时。等值线图如下所示：T=0时（初始时刻）：T=7.6h时：T=15.2h时：计算热物理大作业10T=22.8h时：T=30.4h时：T=38.0h时（稳态时刻）：计算热物理大作业11从上面的图可以看出：（1）刚开始，同水平线温度相同，同一竖直线上温度呈线性分布，符合初始条件的要求。（2）在时间间隔相同的情况下，刚刚开始时温度变化比较迅速，随着时间增加，温度变化趋向缓慢，最后达到稳定。（3）由于左右边界条件对称，理论上温度分布也是左