注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

AlexYoung

做好大家都必须要做的事,发展自己感兴趣的事,结束圣人所不齿的事...

 
 
 

日志

 
 

三维圆柱芯块的温度场数值计算  

2011-04-05 20:24:09|  分类: 课程练习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

%              圆柱芯块的温度场数值计算(堆热工流体数值计算课程报告)
%                                     ——杨万奎     2010.12.10加速
%S.I.P三维柱坐标系,Peric书上的符号系统,初始为恒温进行迭代;相关说明如下:
%分别设定三个空间变量的离散度,依次为半径、高度和角度。我的差分方程建立在了
%除中轴线的所有节点上,因此会出现一些边界节点为已知温度的情况,此时将主对角
%元设置为1,其余对角元设置为零,q(l)设置为对应的环境温度。参数设置如下:
%                    R=0.0045 m-----------半径长度;
%                    Z=0.01 m-------------芯块总高;
%                    K=2.3 W/(m*K)------芯块热导率;
%                    q_v=4*10^8 W/m^3---体积释热率;
%                    T0=400 摄氏度 -------环境温度;
%符号说明如下:
%      P_j-----------绘图时选择的Z=P_j*delta_z位置处的平面的温度分布;
%      P_k-----------绘图时选择的phi=P_k*delta_phi位置处的平面的温度分布;
clear;clc;tic
NI=10;NJ=10;NK=10;
NIJK=NI*NJ*NK;
delta_r=0.0045/NI;delta_z=0.01/(NJ-1);delta_phi=2*pi/NK;
alpha=0.96;
P_j=NJ/2;
P_k=NI/2;
%==========================================================================
%给变量预分配存储空间之后速度会增加非常多,节点越多这种优势越明显
%比如节点取为50*50*50:
%没有以下语句时所需的时间为:1308 seconds
%加上以下语句时所需的时间为:39.674673 seconds.
%    A10~A0~B10---系数矩阵的非零对角元(已知);
%    a10,a6,a3,a1,a0----------L阵的非零对角元;
%    b10,b6,b3,b1-------------U阵的非零对角元;
%    T----------------------温度变量(待求量);
%    q--------------------方程右端值(已知量);
%    R=q-AT------------------------------残差;
%    D=T^(n+1)-T^n-----同一节点两次计算的温度差;
%    y=UD-----------------------------中间变量;
%    X1,Y1,T1;X2,Y2,T2;为绘图用到的坐标及函数值;
A10=zeros(NIJK,1);A6=zeros(NIJK,1);A3=zeros(NIJK,1);A1=zeros(NIJK,1);
B10=zeros(NIJK,1);B6=zeros(NIJK,1);B3=zeros(NIJK,1);B1=zeros(NIJK,1);
A0=zeros(NIJK,1);a0=zeros(NIJK,1);
a10=zeros(NIJK,1);a6=zeros(NIJK,1);a3=zeros(NIJK,1);a1=zeros(NIJK,1);
b10=zeros(NIJK,1);b6=zeros(NIJK,1);b3=zeros(NIJK,1);b1=zeros(NIJK,1);
T=zeros(NIJK,1);
q=zeros(NIJK,1);
R=zeros(NIJK,1);
D=zeros(NIJK,1);
y=zeros(NIJK,1);
X1=zeros(NK+1,NI+1);Y1=zeros(NK+1,NI+1);T1=zeros(NK+1,NI+1);
X2=zeros(NJ,NI+1);  Y2=zeros(NJ,NI+1);  T2=zeros(NJ,NI+1);
%==========================================================================
%初始化,分步进行初始化
%1.将所有对角线扩展为NIJK维的,并初始化为零;
for k=1:NK
    for j=1:NJ
        for i=1:NI
            l=(k-1)*NI*NJ+(j-1)*NI+i;
            r(i)=i*delta_r;
            A10(l)=0;
            A6(l)=0;
            A3(l)=0;
            A1(l)=0;
            A0(l)=1/(r(i)*delta_r)-2/delta_r^2-2/delta_z^2-2/(r(i)*delta_phi)^2;
            B1(l)=0;
            B3(l)=0;
            B6(l)=0;
            B10(l)=0;
        end
    end
end
%2.分别对各条对角线元素进行初始化,赋初值;
for k=1:NK
    for j=1:NJ
        for i=1:NI
            l=(k-1)*NI*NJ+(j-1)*NI+i;
            q(l)=-4e8/2.3;
            if l>1
                A1(l)=1/delta_r^2-1/(r(i)*delta_r);
            end
            if l>NI
                A3(l)=1/delta_z^2;
            end
            if l>NI*NJ
                A6(l)=1/(r(i)*delta_phi)^2;
            end
            if l>NI*NJ*(NK-1)
                A10(l)=1/(r(i)*delta_phi)^2;
            end
            if l+1<NIJK+1
                B1(l)=1/delta_r^2;
            end
            if l+NI<NIJK+1
                B3(l)=1/delta_z^2;
            end
            if l+NI*NJ<NIJK+1
                B6(l)=1/(r(i)*delta_phi)^2;
            end
            if l+NI*NJ*(NK-1)<NIJK+1
                B10(l)=1/(r(i)*delta_phi)^2;
            end
           
            if  (i==NI|j==1|j==NJ)
                A0(l)=1;
                A1(l)=0;
                A3(l)=0;
                A6(l)=0;
                A10(l)=0;
                B1(l)=0;
                B3(l)=0;
                B6(l)=0;
                B10(l)=0;
                q(l)=400;
            end
        end
    end
end
%============================================
%求L、U矩阵元并同时可以解下三角阵
%迭代过程中的符号为LU*D=R,其中:
%    R=q-AT       -----------------------残差;
%    D=T^(n+1)-T^n-----同一节点两次计算的温度差;
%    y=UD;
for   l=1 : NIJK
        if l>1
            b11=b1(l-1);
            b31=b3(l-1);
            b61=b6(l-1);
            b101=b10(l-1);
        else
            b11=0;
            b31=0;
            b61=0;
            b101=0;
        end
        if l>NI
            b1NI=b1(l-NI);
            b3NI=b3(l-NI);
            b6NI=b6(l-NI);
            b10NI=b10(l-NI);
        else
            b1NI=0;
            b3NI=0;
            b6NI=0;
            b10NI=0;
        end
        if l>NI*NJ
            b1NINJ=b1(l-NI*NJ);
            b3NINJ=b3(l-NI*NJ);
            b6NINJ=b6(l-NI*NJ);
            b10NINJ=b10(l-NI*NJ);
        else
            b1NINJ=0;
            b3NINJ=0;
            b6NINJ=0;
            b10NINJ=0;
        end
        if l>NI*NJ*(NK-1)
            b1NINJNKM1=b1(l-NI*NJ*(NK-1));
            b3NINJNKM1=b3(l-NI*NJ*(NK-1));
            b6NINJNKM1=b6(l-NI*NJ*(NK-1));
            b10NINJNKM1=b10(l-NI*NJ*(NK-1));
        else
            b1NINJNKM1=0;
            b3NINJNKM1=0;
            b6NINJNKM1=0;
            b10NINJNKM1=0;
        end
          a10(l)=A10(l)/(1+alpha*(b1NINJNKM1+b3NINJNKM1+2*b6NINJNKM1));
          a6(l)=A6(l)/(1+alpha*(b1NINJ+b3NINJ));
          a3(l)=A3(l)/(1+alpha*(b1NI+b6NI+b10NI));
          a1(l)=A1(l)/(1+alpha*(b31+b61+b101));
          a0(l)=A0(l)-(a1(l)*b11+a3(l)*b3NI+a6(l)*b6NINJ+a10(l)*b10NINJNKM1)+alpha*...
               (a10(l)*b1NINJNKM1+a10(l)*b3NINJNKM1+a10(l)*b6NINJNKM1+...
                a6(l)*b1NINJ+a6(l)*b3NINJ+a3(l)*b1NI+...
                a1(l)*b31+a1(l)*b61+a3(l)*b6NI+...
                a1(l)*b101+a3(l)*b10NI+a6(l)*b10NINJ);
          b1(l)=(B1(l)-alpha*(a3(l)*b1NI+a6(l)*b1NINJ+a10(l)*b1NINJNKM1))/a0(l);
          b3(l)=(B3(l)-alpha*(a1(l)*b31+a6(l)*b3NINJ+a10(l)*b3NINJNKM1))/a0(l);
          b6(l)=(B6(l)-alpha*(a1(l)*b61+a3(l)*b6NI))/a0(l);
          b10(l)=(B10(l)-alpha*(a1(l)*b101+a3(l)*b10NI+2*a6(l)*b10NINJ))/a0(l);
          T(l)=1000;
end
%==========================================================================
%进行内迭代
%收敛准则:delta<=10,为防止发散时无谓的计算,取delta>400就停止迭代。
delta=1;
n=0;
while (delta>1e-5)
    %---------------------------------------------
    %求残差,同时解下三角阵的内迭代
    for  l=1 : NIJK
            if  l-1>0
                TM1=T(l-1);
                y1=y(l-1);
            else
                TM1=0;
                y1=0;
            end
            if  l-NI>0
                TMNI=T(l-NI);
                yNI=y(l-NI);
            else
                TMNI=0;
                yNI=0;
            end
            if  l-NI*NJ>0
                TMNINJ=T(l-NI*NJ);
                yNINJ=y(l-NI*NJ);
            else
                TMNINJ=0;
                yNINJ=0;
            end
            if  l-NI*NJ*(NK-1)>0
                TMNINJNKM1=T(l-NI*NJ*(NK-1));
                yNINJNKM1=y(l-NI*NJ*(NK-1));
            else
                TMNINJNKM1=0;
                yNINJNKM1=0;
            end
            %----------------
            if  l+1<NIJK+1
                TP1=T(l+1);
            else TP1=0;
            end
            if  l+NI<NIJK+1
                TPNI=T(l+NI);
            else TPNI=0;
            end
            if  l+NI*NJ<NIJK+1
                TPNINJ=T(l+NI*NJ);
            else TPNINJ=0;
            end
            if  l+NI*NJ*(NK-1)<NIJK+1
                TPNINJNKM1=T(l+NI*NJ*(NK-1));
            else TPNINJNKM1=0;
            end
            R(l)=q(l)-(A10(l)*TMNINJNKM1+A6(l)*TMNINJ+A3(l)*TMNI+A1(l)*TM1+A0(l)*T(l)...
                +B1(l)*TP1+B3(l)*TPNI+B6(l)*TPNINJ+B10(l)*TPNINJNKM1);
            y(l)=(R(l)-(a1(l)*y1+a3(l)*yNI+a6(l)*yNINJ+a10(l)*yNINJNKM1))/a0(l);
    end
    %-----------------------------------------------------
    %解上三角阵的内迭
    for l=NIJK:-1:1
           if  l+1<NIJK+1
               y1=D(l+1);
           else
               y1=0;
           end
           if  l+NI<NIJK+1
               yNI=D(l+NI);
           else
               yNI=0;
           end
           if  l+NI*NJ<NIJK+1
               yNINJ=D(l+NI*NJ);
           else
               yNINJ=0;
           end
           if  l+NI*NJ*(NK-1)<NIJK+1
               yNINJNKM1=D(l+NI*NJ*(NK-1));
           else
               yNINJNKM1=0;
           end
           D(l)=y(l)-(b1(l)*y1+b3(l)*yNI+b6(l)*yNINJ+b10(l)*yNINJNKM1);
           T(l)=T(l)+D(l);
           d(l)=D(l)/T(l);
    end
    delta=max(abs(d))
    n=n+1;
end
%==========================================================================
Tnmax=max(T)
Tnmin=min(T)
n
delta
%==========================================================================
%绘图,取z为某常数的平面绘制三维温度分布绘图
s=0;
for k=1:NK
    for i=2:NI+1
        l=(k-1)*NI*NJ+(P_j-1)*NI+i-1;
        X1(k,i)=(i-1)*delta_r*cos(k*delta_phi);
        Y1(k,i)=(i-1)*delta_r*sin(k*delta_phi);
        T1(k,i)=T(l);
    end
    s=s+2*T1(k,2)-T1(k,3);
end
Tzmax=max(T1(:))
[a,b]=find(T1==Tzmax)
%处理中心点的绘图,取为周围点的平均值;
for k=1:NK+1
    X1(k,1)=0;
    Y1(k,1)=0;
    T1(k,1)=s/NK;
end
%处理缺口的绘图,因为绘图时不会自动封闭,故需要将角度为0的值赋值给2Pi处;
for i=2:NI+1
    X1(NK+1,i)=(i-1)*delta_r*cos(1*delta_phi);
    Y1(NK+1,i)=(i-1)*delta_r*sin(1*delta_phi);
    T1(NK+1,i)=T1(1,i);
end
surf(X1,Y1,T1)
%view(-90,0)
%--------------------------------------------------------------------------
%绘制r-z平面的温度分布绘图
figure
for j=1:NJ
    for i=2:NI+1
        l=(P_k-1)*NI*NJ+(j-1)*NI+i-1;
        X2(j,i)=(i-1)*delta_r;
        Y2(j,i)=(j-1)*delta_z;
        T2(j,i)=T(l);
    end
    X2(j,1)=0;
    Y2(j,1)=(j-1)*delta_z;
    T2(j,1)=2*T2(j,2)-T2(j,3);
end
surf(X2,Y2,T2)
%==========================================================================
toc

 

绘图结果及分析:

设置项:网格划分:50×50×50;

收敛标准:delta<1e-5

Elapsed time is 37.375709 seconds.

绘图控制命令:

view(0,90)
axis([-0.0045,0.0045,-0.0045,0.0045])
colorbar
grid off

三维圆柱芯块的温度场数值计算 - 青楼薄名 - Alex小奎子
 

 

绘图控制命令:

view(0,90)
axis([0,0.0045,0,0.01])
colorbar
grid off

三维圆柱芯块的温度场数值计算 - 青楼薄名 - Alex小奎子

绘图控制命令:

shading interp

三维圆柱芯块的温度场数值计算 - 青楼薄名 - Alex小奎子

 

三维圆柱芯块的温度场数值计算 - 青楼薄名 - Alex小奎子

 

  评论这张
 
阅读(476)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018