《计算热物理(共30页).docx》由会员分享，可在线阅读，更多相关《计算热物理(共30页).docx（30页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上目录244569121416题目：定义在正方形区域（如图所示）H*H=3m*3m 的二维非稳态导热问题的控制方程及其定解条件为其中，c，分别为导热体的密度、比热和导热系数，且=7820kg/m3，c=460J/(kgK)，=15W/(mK)。将定义域各方向均匀分成10个区块，取点中心网格分割方式，各边界共 11个节点; 或者取块中心网格分割，各边界共 10个节点，（两种网格，至少做一种），按以下要求，数值求解该问题。 （1） 取（X，Y）=（x，y）H，=T-T0T1-T0，=tcH2/，将方程和定解条件无量纲化，并在无量纲化条件下求解。（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的交替方向线迭代。并对收敛结果做出比
3、较。（5） 选择全隐式格式在 5个不同时刻（包括 1个到达稳态的时刻）的计算结果，将无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。从 5个不同时刻的等值线图，分析时间发展中 T的变化特性。（6）根据数值计算情况，（精度高低，CPU时间消耗，编制程序难易等）对各种不同的离散方法及代数解法进行分析比较，并做出结论。（7）综合上述内容和作业情况，写出作业报告，并对此类作业提出建议和要求。解：（1）方程和定解条件的无量纲化：将（X，Y）=（x，y）H，=T-T0T1-T0，=tcH2/代入原方程组，无量纲化，得到：（2） 控制方程的离散：A 控制容积法假定非稳态项中温度随空间为阶梯分
4、布，扩散项中温度随时间线性分布、随空间阶梯分布显式格式: P-P0XY=E0-P0Xe-P0-W0XwY+N0-P0Yn-P0-S0YsX 写成：aPP=aEE0+aWW0+aNN0+aSS0+b 的形式， 则：aE=YXe，aW=YXw，aN=XYn，aS=XYs aP=XY，aP0=XY-aE-aW-aN-aS，b=aP0P0 全隐式格式： P-P0XY=E-PXe-P-WXwY+N-PYn-P-SYsX 写成：aPP=aEE+aWW+aNN+aSS+b 的形式， 则：aE=YXe，aW=YXw，aN=XYn，aS=XYs aP=XY+aE+aW+aN+aS，aP0=XY，b=aP0P0B
5、 ADI格式：i,j*-i,jn/2=i+1,j*-2i,j*+i-1,j*X2+i,j+1n-2i,jn+i,j-1nY2 i,jn+1-i,j*/2=i+1,j*-2i,j*+i-1,j*X2+i,j+1n+1-2i,jn+1+i,j-1n+1Y2（3）稳定性分析:A 显示格式：aPP=aEE0+aWW0+aNN0+aSS0+baPg=aEeikXX+aWe-ikXX+aNeikYY+aSe-ikYY+aP0 (*)本题中定义域各方向均匀分成10个区块，则Xe=Xw=Yn=Ys=X=Y=0.1aE=aW=aN=aS=1，aP=XY则(*)式可化为：g=2coskXX+2coskYY+XY/
6、-4XY/假设kX=kY，则g=1-8X2sin2kX2有g1恒成立，则X2sin2kX214X214，14X2。即14X2时，该显式格式收敛。本题中，取X=0.1，故0.0025。B 全隐式格式：aPP=aEE+aWW+aNN+aSS+baPg=gaEeikXX+aWe-ikXX+aNeikYY+aSe-ikYY+aP0(*)本题中定义域各方向均匀分成10个区块，则Xe=Xw=Yn=Ys=X=Y=0.1aE=aW=aN=aS=1，aP0=XY(*)式化为：g=XY/4-2coskXX-2coskYY+XY/显然g1恒成立，故全隐式格式恒稳定。C ADI格式：i,j*-i,jn/2=i+1,j
7、*-2i,j*+i-1,j*X2+i,j+1n-2i,jn+i,j-1nY2 i,jn+1-i,j*/2=i+1,j*-2i,j*+i-1,j*X2+i,j+1n+1-2i,jn+1+i,j-1n+1Y2空间导数离散式用微分算子符号写出，上面两个等式可化为：1-X2X2i,j*=1+Y2Y2i,jn，X=/X21-Y2Y2i,jn+1=1+X2X2i,j*，Y=/Y2消去i,j*，得：1-X2X21-Y2Y2i,jn+1=1+X2X21+Y2Y2i,jn故g=1-2Xsin2kXX21-2Ysin2kYY21+2Xsin2kXX21+2Ysin2kYY2显然g1恒成立，故ADI格式恒稳定。（4
8、） 温度场的计算求解将定义域各方向均匀分成 10 个区块,采用点中心网格分割方式，各边界共11 个节点。边界条件的离散：（1）Y=0：i,1=1(i=111)；（2）Y=1：i,11=0(i=111)；（3）X=0：j=210显式：a1,j1,j=a2,j2,j+a1,j+11,j+1+a1,j-11,j-1+b其中： a1,j=XY2，a2,j=YXe，a1,j+1=X2Yn，a1,j-1=X2Ys b=1,j0a1,j-a2,j-a1,j+1-a1,j-1+qwY隐式：a1,j1,j=a2,j2,j+a1,j+11,j+1+a1,j-11,j-1+b其中： a2,j=YXe，a1,j+1=
9、X2Yn，a1,j-1=X2Ys a1,j=XY2+ a2,j+a1,j+1+a1,j-1，b=XY21,j0+qwY（4）X=1：j=210显式：a11,j11,j=a10,j10,j+a11,j+111,j+1+a11,j-111,j-1+b其中：a11,j=XY2，a10,j=YXw，a11,j+1=X2Yn，a11,j-1=X2Ysb=11,j0a11,j-a10,j-a11,j+1-a11,j-1+qeY隐式：a11,j11,j=a10,j10,j+a11,j+111,j+1+a11,j-111,j-1+b其中： a10,j=YXw，a11,j+1=X2Yn，a11,j-1=X2Ys
10、 a11,j=XY2+ a10,j+a11,j+1+a11,j-1，b=XY211,j0+qeY初始条件：由“初始T(x,y,0)为x=0x=H的温度成线性分布”我们可以推知：i,j=1-0.1*(i-1) i,j=111据此编程，编程软件：Visual C+，程序代码见附录。其中，左右边界的热流量qe和qw、时间步长、稳态精度、迭代精度、时间上限以及最大迭代次数均可在define项目中直接调整。几种迭代算法的比较：1. Jacobi简单点迭代：利用上一时刻的节点温度值计算下一时刻该节点的温度，程序较为简单。收敛速度慢，迭代最大次数为9次，出现在时间刚开始的阶段。下面是计算结果：2. Gaus
11、s-Seidel点迭代：在计算节点值时，采用邻近节点的最新可用值，而不局限于前一迭代完成后所得到的值。此种方法由于及时引入新值，迭代收敛速度快于简单迭代。在实际程序设计中，先计算边界，使边界的最新可用值被尽可能快的引入到内节点的计算中，对加快收敛速度起到很好的作用。在以后的程序中，也都采用了这种方法。该算法的最大迭代次数为7次，出现在时间值很小的起步阶段。下面是计算结果：3. 松弛迭代本次上机取松弛因子w在从0.1开始按 0.1增量递增到1.9。松弛迭代的收敛速度与松弛因子有关。当松弛因子在1.2-1.4时，收敛速度明显加快，而当松弛因子取0.1或1.9附近的值时，有可能发散。4. Gauss
12、-Seidel的交替方向线迭代此法扫描方向可以变化，且有多种组合。本程序先处理边界条件，然后采用从左至右的列扫描，后再采用从上至下的行扫描，全场两次扫描组成一轮迭代，这种扫描是对一条线上隐式格式迭代求解。相对于Gauss-Seidel迭代，收敛速度有所提高，最大迭代次数为4次。下面是计算结果：（5）温度场的等值线图这里，结合 Gauss-Seidel的交替方向线迭代的计算结果，无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。我们可以看到，达到稳态所需的无量纲时间为0.636，折算成实际时间为38小时。等值线图如下所示：T=0时（初始时刻）：T=7.6h时：T=15.2h时：T=
13、22.8h时：T=30.4h时：T=38.0h时（稳态时刻）：从上面的图可以看出： （1）刚开始，同水平线温度相同，同一竖直线上温度呈线性分布，符合初始条件的要求。 （2）在时间间隔相同的情况下，刚刚开始时温度变化比较迅速，随着时间增加，温度变化趋向缓慢，最后达到稳定。 （3）由于左右边界条件对称，理论上温度分布也是左右对称。数值解得到的结果也验证了这一点。最高温度为1.0885，折算成实际温度为413.275，出现在图的左下角和右下角。（4）随着时间的推进，低温区逐渐缩小，高温区逐渐扩大。与左右边界有恒定热流进入此区域有关。 综上可见，计算结果与实际情况相符。（6）分析比较这里先把显示格式和
14、ADI格式的计算结果予以展示，以作为对前面隐式格式迭代求解计算结果的对比。显示格式计算结果：ADI格式计算结果：A 精度从几种算法的运算结果来看，精度上的差别不大；但从达到稳态所需时间来看，精度上还是有差别。计算结果的精度主要取决于离散方程的精度，由离散方程截断误差的量级所决定。本题中所使用的各种迭代算法的精度均为二阶精度。时间方面，由于求解下一时刻温度场时，显示格式和ADI格式上一时刻的温度场并没有收敛，而隐式格式上一时刻的温度场是收敛的，两者的差异导致最终达到稳态所需的时间不同，相差0.07，约为33min。B CPU时间消耗显示格式和ADI格式的求解不需要迭代，运算量小，所以耗时短，分别
15、为0.062秒和0.078秒。相比之下，隐式格式的算法普遍耗时长。耗时最多的是Gauss-Seidel迭代，这是因为相比于Jacobi迭代，前者对边界的处理采用了一种更容易理解的方式，但需要额外的空间，也增加了运算量，所以虽然迭代次数有所减少，但总的运行时间反而增加。Gauss-Seidel线迭代的运行时间明显比Gauss-Seidel迭代减少很多，这主要归功于迭代次数的进一步减少。C 编程难易程度显示格式的程序设计起来简单，不需要考虑迭代等一系列问题。ADI格式也是直接求解，程序较为简单。隐式格式的求解，相对前两者就复杂了一些，各种算法间也有较大差异。Jacobi迭代和Gauss-Seide
16、l迭代算法都较简单，编程难度几乎同等，Gauss-Seidel线迭代的算法相对难一些，编程难度也大了一些。而带松弛因子的迭代，需要在时间步长和迭代次数之上再增加一个变量，程序设计的难度大大增加，其中带松弛的Gauss-Seidel线迭代的程序部分是本人认为最难设计的。 结论：显示和ADI格式，虽然可以直接求解，但精度不高，受时间步长的影响也很大。带松弛的迭代需要考虑松弛因子的选取，若选取不当将会带来巨大的运算量，而且程序设计难度大。隐式格式不受时间步长的影响，Gauss-Seidel线迭代的迭代次数少，CPU运行时间短，程序设计难度可以接受。综合考虑，本人认为隐式格式Gauss-Seidel线
17、迭代是最为理想的求解方法。（7）个人心得体会此次计算热物理大作业的完成，前后用了将近一周的时间，难度最大的部分在于程序源代码的编写。在代码的编写工程中，我们需要考虑多种参量与变量，以及各种各样的算法。其中，最难处理的部分在于边界条件，这不仅仅因为需要考虑的因素较内节点多，还因为需要编写专门的算法，因而增加了程序的设计难度，而且，稍有不慎就满盘皆输。通过好几百行代码的编写，巩固了既有的程序设计基础，又学到了一些新的东西，如程序执行时间CPU耗时的计算。就计算热物理这门学科而言，通过这次大作业，将课堂所学加以实际应用，学以致用，不仅扎实了基础，还得到了一些感悟和收获。最后，说一点个人的看法：二维数
18、组的求解，最起码的时间复杂度为N的平方项，如果考虑迭代，增加迭代次数，则又增加了一阶复杂性，现在我们所处理的非稳态问题，本身就需要考虑时间的递增，至此，我们可以保守估计，程序的时间复杂度已经达到了N的4次方。如果再加入可变的松弛因子，复杂度就达到5次方。时间复杂度，其实也代表了程序设计的复杂性，一个5重循环就足以让人头疼不已，更何况还得考虑多种复杂变化的变量。所以，我建议在大作业中舍去松弛因子这一项，使复杂度降阶。如果要让我们认识到松弛因子的作用以及不同松弛因子对迭代次数的影响，可以专门设计一个不包含时间项的问题。附录：（源代码部分）1. 显式格式：#include#include#inclu
19、de#define dx 0.1 /x方向步长#define dy 0.1 /y方向步长#define qw 1.0#define qe 1.0#define n (int)(1/dx)+2 /n=12#define tup 0.002 /时间步长，不大于0.0025#define tup_limit 50 /时间上限#define epsilon 0.00001using namespace std;/数组初始化void init(double Tnn)int i,j;for(i=1;in;i+)for(j=1;j0;i-)for(j=1;jn;j+)coutsetw(5)Tij ;cout
20、endl;/绝对值函数double re(double x)if(x0.0)x=0.0-x;return x;/向量范数int fanshu(double X,double Y)int i;double s=0.0;for(i=1;in;i+)if(sre(Xi-Yi)s=re(Xi-Yi);if(sepsilon)return 0;else return 1;/计算温度值void fun2(double Tnn,double tempnn)int i,j;double c1,c2;c1=dx*dy/tup;for(i=2;in-1;i+) /内部各点 i,j=2-10for(j=2;jn-1
21、;j+)c2=tempij-1+tempij+1+tempi-1j+tempi+1j;c2+=tempij*(c1-4.0);Tij=c2/c1;/左边界,TDMA方法求解c1=c1/2.0;temp01=0.0;for(i=2;in-1;i+)temp0i=0.5/(c1-0.5*temp0i-1);temp12=(temp21*(c1-2)+qw*dy+T22+0.5*T11)/c1;for(i=3;i1;i-)Ti1=temp0i*Ti+11+temp1i;/右边界,TDMA方法求解for(i=2;in-1;i+)temp0i=0.5/(c1-0.5*temp0i-1);temp12=(
22、temp211*(c1-2)+qe*dy+T210+0.5*T111)/c1;for(i=3;i1;i-)Ti11=temp0i*Ti+111+temp1i;/时间向前走一步，并判断是否稳态/与前一时刻的温度进行比较，若低于精度值则返回0，否则返回1int fun1(double Tnn)double tempnn;int i,j;for(i=1;in;i+)for(j=1;jn;j+)tempij=Tij; /保存前一时刻的值fun2(T,temp); /计算该时刻的值for(i=2;in-1;i+) /逐行比较温度值if(fanshu(Ti,tempi)!=0)return 1;retur
23、n 0;/主函数int main()clock_t start,end;start = clock(); /程序运行时间double Tnn;init(T); /数组初始化/print(T); /打印数组double time=0.0; /无量纲时间while(time=tup_limit)time+=tup;if(fun1(T)=0)cout无量纲时间time时达到稳态endl;cout稳态时区域各节点的无量纲温度为:endl;coutsetiosflags(ios:fixed)setprecision(4);print(T);end = clock();cout运行时间:(double)(
24、end - start) / CLOCKS_PER_SECSendl;return 0; /结束程序运行cout在无量纲时间tup_limit内没有达到稳态.endl;end = clock();cout运行时间(double)(end - start) / CLOCKS_PER_SEC秒endl;return 0;2. ADI格式：#include#include#include#define dx 0.1 /x方向步长#define dy 0.1 /y方向步长#define qw 1.0#define qe 1.0#define n (int)(1/dx)+2 /n=12#define t
25、up 0.002 /时间步长，不大于0.0025#define tup_limit 10 /时间上限#define epsilon 0.00001using namespace std;/数组初始化void init(double Tnn)int i,j;for(i=1;in;i+)for(j=1;j0;i-)for(j=1;jn;j+)coutsetw(5)Tij ;coutendl;/绝对值函数double re(double x)if(x0.0)x=0.0-x;return x;/向量范数int fanshu(double X,double Y)int i;double s=0.0;fo
26、r(i=1;in;i+)if(sre(Xi-Yi)s=re(Xi-Yi);if(sepsilon)return 0;else return 1;/计算温度值void fun2(double Tnn,double tempnn)int i,j;double c1;double p9,q9;/左边界,TDMA方法求解c1=dx*dy/tup;p0=1.0/c1;for(i=1;i8;i+)pi=1.0/(c1-pi-1);q0=(T22*2.0+2.0*qw*dy+T11+(c1-4)*temp21)/c1;for(i=1;i1;i-)Ti1=pi-2*Ti+11+qi-2; /右边界,TDMA方
27、法求解q0=(T210*2.0+2.0*qe*dy+T111+(c1-4)*temp211)/c1;for(i=1;i1;i-)Ti11=pi-2*Ti+111+qi-2;/内部各点ADI算法c1=2.0*dx*dy/tup;p0=1.0/(c1+2.0);for(i=1;i8;i+)pi=1.0/(c1+2.0-pi-1);for(j=2;j11;j+) /列扫描 q0=(T2j*(c1-2)+T2j+1+T2j-1+T1j)/(c1+2.0); for(i=1;i1;i-)Tij=pi-2*Ti+1j+qi-2;/行扫描for(i=2;i11;i+) q0=(Ti2*(c1-2)+Ti+1
28、2+Ti-12+Ti1)/(c1+2.0); for(j=1;j1;j-)Tij=pj-2*Tij+1+qj-2;/时间向前走一步，并判断是否稳态/与前一时刻的温度进行比较，若低于精度值则返回0，否则返回1int fun1(double Tnn)double tempnn;int i,j;for(i=1;in;i+)for(j=1;jn;j+)tempij=Tij; /保存前一时刻的值fun2(T,temp); /计算该时刻的值for(i=2;in-1;i+) /逐行比较温度值if(fanshu(Ti,tempi)!=0)return 1;return 0;/主函数int main()cloc
29、k_t start,end;start = clock(); /程序运行时间double Tnn;init(T); /数组初始化print(T); /打印数组double time=0.0; /无量纲时间while(time=tup_limit)time+=tup;if(fun1(T)=0)cout无量纲时间time时达到稳态endl;cout稳态时区域各节点的无量纲温度为:endl;coutsetiosflags(ios:fixed)setprecision(4);print(T);end = clock();cout运行时间:(double)(end - start) / CLOCKS_P
30、ER_SECSendl;return 0; /结束程序运行cout在无量纲时间tup_limit内没有达到稳态.endl;end = clock();cout运行时间(double)(end - start) / CLOCKS_PER_SEC秒endl;return 0;3. Jacobi简单点迭代：#include#include#include#define dx 0.1 /x方向步长#define dy 0.1 /y方向步长#define qw 1.0#define qe 1.0#define n (int)(1/dx)+2 /n=12#define tup 0.002 /时间步长，不大
31、于0.0025#define tup_limit 10 /时间上限#define epsilon 0.00001 /稳态精度#define epsilon2 0.00001 /迭代精度#define nmax 10000 /最大迭代次数using namespace std;/数组初始化void init(double Tnn)int i,j;for(i=1;in;i+)for(j=1;j0;i-)for(j=1;jn;j+)coutsetw(5)Tij ;coutendl;/绝对值函数double re(double x)if(x0.0)x=0.0-x;return x;/向量范数int f
32、anshu(double X,double Y)int i;double s=0.0;for(i=1;in;i+)if(sre(Xi-Yi)s=re(Xi-Yi);if(sepsilon)return 0;else return 1;int Jacobi(double Tnn,double Tmidnn,double tempnn)int i,j;double c1=dx*dy/tup,c2;/左、右边界for(j=2;jn-1;j+)c2=2.0*Tmidj2+Tmidj+11+Tmidj-11+c1*tempj1+2.0*qw*dy; Tj1=c2/(c1+4.0); c2=2.0*Tmi
33、dj10+Tmidj+111+Tmidj-111+c1*tempj11+2.0*qe*dy; Tj11=c2/(c1+4.0);/内部节点for(i=2;in-1;i+)for(j=2;jn-1;j+)c2=Tmidij-1+Tmidij+1+Tmidi+1j+Tmidi-1j+c1*tempij;Tij=c2/(c1+4.0);/判断是否收敛for(i=2;in-1;i+)for(j=1;jepsilon2)return 0;return 1;/计算温度值int fun2(double Tnn,double tempnn)int i,j;int m=0;double Tmidnn;while
34、(m=nmax)for(i=1;in;i+) /保存上一次迭代后的值for(j=1;jn;j+)Tmidij=Tij;if(Jacobi(T,Tmid,temp)!=0)return m+1;else m+=1;return m+1;/时间向前走一步，并判断是否稳态/与前一时刻的温度进行比较，若低于精度值则返回0，否则返回1int fun1(double time,double Tnn)double tempnn;int i,j,m;for(i=1;in;i+)for(j=1;jnmax)cout时刻time时矩阵不收敛，以下计算不正确!endl;else if(int(time*500)%10=0)cout求解时刻time时温度场所用迭代次数为:mendl;for(i=2;in-1;i+) /逐行比较温度值if(fanshu(Ti,tempi)!=0)return 1;return 0;/主函数int main()clock_t start,end;start = clock(); /程序运行时间double Tnn;init(T); /数组初始化print(T); /打印数组double time=0.0; /无量纲时间while(time=tup_limit)time+=tup;if(fun1(time,T)=0)cout无量纲时间time时达到稳态endl;