数值计算求解静态热传导方程

首先要用到一个三维画图的函数，在三维空间中用三角形连接。

function quadplot(nodes,elements,sol)
Zeile=size(elements,1);
x=nodes(:,1);
y=nodes(:,2);
z=sol;
T=zeros(2*Zeile,3);
%把平面上的一个四边形elements切割成2个三角形，所以是2*Zeile,再转化到三维空间去,一个四边形的elements包含四个点，切割成两个三角形之后，每个三角形包含三个点。
 for i=1:size(elements,1)
    T(2*i-1,:)=[elements(i,1:2),elements(i,4)];
    T(2*i,:)=[elements(i,2:4)];
 end
trisurf(T,x,y,z);
xlabel('X'),ylabel('Y'),zlabel('Z');
end
%elements = [1, 2, 5, 4;
             2, 3, 6, 5; 
             4, 5, 8, 7; 
             5, 6, 9, 8]
%T=[1 2 5
    5 4 1
    2 3 6
    6 5 2
    4 5 8
    4 8 7
    5 6 9
    5 8 9]

先看看静态热传导方程的函数表达式，看第七章的课件第15页，可以看到这里要用到之前的linquadderivref(xi,eta), getJacobian(nodes,xi,eta), gx2dref, gw2dref,来复习一下这几个函数

function deriv=linquadderivref(xi,eta)
         deriv=[-1/4*(1-eta),-1/4*(1-xi)
                1/4*(1-eta), -1/4*(1+xi)
                1/4*(1+eta),1/4*(1+xi)
                -1/4*(1+eta),1/4*(1-xi)];
end

function[J,detJ,invJ]=getJacobian(nodes,xi,eta)
Jt=linquadderivref(xi,eta)'*nodes;
J=Jt';
detJ=det(J);
invJ=inv(J);
end

function gaussx = gx2dref(n)
if n==1
   gaussx=[0 0];
elseif n==2
   gaussx=[-1/sqrt(3) -1/sqrt(3)
           -1/sqrt(3) 1/sqrt(3)
           1/sqrt(3)  -1/sqrt(3)
           1/sqrt(3)  1/sqrt(3)];
elseif n==3
   gaussx=[-sqrt(3/5) -sqrt(3/5) 
            -sqrt(3/5) 0 
            -sqrt(3/5) sqrt(3/5)
            0          -sqrt(3/5) 
            0          0 
            0          sqrt(3/5) 
            sqrt(3/5)  -sqrt(3/5)  
            sqrt(3/5)  0
            sqrt(3/5)  sqrt(3/5) ];
else
    disp('false n')
end

function gaussw = gw2dref(n)
if n==1
    gaussw=[4];
elseif n==2
    gaussw=[1
            1
            1
            1];
elseif n==3
    gaussw=[5/9*5/9 
            5/9*8/9 
            5/9*5/9
            8/9*5/9 
            8/9*8/9
            8/9*5/9
            5/9*5/9 
            5/9*8/9 
            5/9*5/9];
else
    disp('false n')
end

看函数表达式，我们先要做的是将elements一个一个求解出来，看课件的讲解，一个element有四个点，所以我们求出来的elemat是一个4x4的矩阵

function [elemat,elevec]=evaluate_stat(elenodes,gx2dref,gw2dref)
elemat=zeros(4,4);
%按照公式右边，求出来的是一个一个的值，然后把一个一个的值代入这个4x4的零矩阵
lamda=48;
%一个常数，题目会给的
for k=1:4
   for i=1:4
       s=0;%求的这些值是一个一个的累加，所以需要设置一个初值
       for j=1:size(gw2dref);%要循环所有的高斯点
           xi=gx2dref(j,1);
           eta=gx2dref(j,2);
           deriv=linquadderivref(xi,eta);
          [J,detJ,invJ]=getJacobian(elenodes,xi,eta);
          Nv=deriv(k,:)*invJ;
          Nt=deriv(i,:)*invJ;
          s=s+lamda*dot(Nv,Nt)*detJ*gpw(j);
      end
   elemat(k,i)=s;
  end
end
elevec=zeros(4,1);
end

下一步要做的就是把这些elemat从lokal变成global，就是进行一下坐标变换

function [sysmat,rhs]=assemble(elemat,elevec,sysmat,rhs,ele]
%ele是local 1234对应到global
a=zeros(size(sysmat));
b=zeros(size(rhs));
for i=1:4
    for j=1:4
    a(ele(i),ele(j))=elemat(i,j);
    end
    b(ele(i))=elevec(i);
end
sysmat=sysmat+a;
rhs=rhs+b
end
%
([1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12; 13, 14, 15, 16], [17; 18; 19; 20], eye(5,5), ones(5,1), [5, 3, 1, 2]) →
[[12,12,10,0,9;15,17,14,0,13;7,8,7,0,5;0,0,0,1,0;3,4,2,0,2],[20;21;19;1;18]]

接着就是处理边界温度了，我们要解的方程是AT=0，解出T矩阵，但如果有边界温度的话，AT就不等于0矩阵了，具体看课件22页举的例子。

function [sysmat,rhs]=assignDBC(sysmat,rhs,dbc)
%dbc是一个两列的矩阵，第一列代表节点位置，第二列代表该节点对应的边界温度
a=size(dbc,1);
%返回该矩阵的行数
for i=1:a
    sysmat(dbc(i,1),:)=0;
    %有边界温度的这一行先全部用零替换
    sysmat(dbc(i,1),dbc(i,1))=1;
    %dbc温度所在的节点用1替换
    rhs(dbc(i,1))=dbc(i,2);
    %把dbc的值放入rhs
end
end
%
([12, 12, 10, 0, 9; 15, 17, 14, 0, 13; 7, 8, 7, 0, 5; 0, 0, 0, 1, 0; 3, 4, 2, 0, 2], [20; 21; 19; 1; 18], [2, −7; 5, −2]) →
[[12,12,10,0,9;0,1,0,0,0;7,8,7,0,5;0,0,0,1,0;0,0,0,0,1],[20;-7;19;1;-2]]

然后就是写这个作业的题目了

先要把所有的节点坐标求出来；然后把所有element的位置求出来；之后就是dbc的边界温度条件。

clear;
r=0.02;
b=0.3;
h=0.3;
nodes=[0 0;b/3 0;2*b/3 0;b 0;
       0 h/3;b/3 h/3;2*b/3 h/3;b h/3;
       0 2*h/3;b/3 2*h/3;2*b/3 2*h/3;b-r*sin(pi/6) h-r*cos(pi/6);
       b h-r;b-r*cos(pi/6) h-r*sin(pi/6);0 h;b/3 h;
       b/2 h;b-r h];
elements=[1 2 6 5;2 3 7 6;3 4 8 7;5 6 10 9;6 7 11 10;11 7 12 14;
          7 8 13 12;9 10 16 15;10 11 17 16;11 14 18 17];
dbc=[1 600;2 600;3 600;4 600;
     12 300;13 300;14 300;18 300];
sysmat=zeros(18,18);
rhs=zeros(18,1);
%一共有18个点组成
for n=1:size(elements,1)
    %一共有size(elements,1)个elements，这里是10个
    [elemat,elevec]=evaluate_stat(nodes(elements(n,:),:),:),gx2dref(2),gw3dref(2));
    %这里是把每个element对应的点代入进行运算
    [sysmat,rhs]=assemble(elemat,elevec,sysmat,rhs,elements(n,:));
    %每个element在local对应的是1 2 3 4个点，但在global就不是了，看给的图就明白了，直接给出了每个element的四个点在global对应的点编号。
end
[sysmat,rhs]=assinDBC(sysmat,rhs,dbc);
T=sysmat\rhs;
%mldivide(sysmat,rhs)
quadplot(nodes,elements,T)
shading interp; grid on
colormap(hot);
colorbar

这个题目的最后一问要让(15, 16, 17, 18）这四个点的温度不超过450K，每次r=r+0.01，可以直接每次r增加0.1代入上面的函数，也可以写一个新的函数进行自动迭代

clear;
r=0.02;b=0.3;h=0.3;T=600*ones(18,1);
while T(15)>450|T(16)>450|T(17)>450|T(18)>450
    nodes=[0 0;b/3 0;2*b/3 0;b 0;
           0 h/3;b/3 h/3;2*b/3 h/3;b h/3;
           0 2*h/3;b/3 2*h/3;2*b/3 2*h/3;b-r*sin(pi/6) h-r*cos(pi/6);
           b h-r;b-r*cos(pi/6) h-r*sin(pi/6);0 h;b/3 h;
           b/2 h;b-r h];
    elements=[1 2 6 5;2 3 7 6;3 4 8 7;5 6 10 9;6 7 11 10;11 7 12 14;
              7 8 13 12;9 10 16 15;10 11 17 16;11 14 18 17];
    dbc=[1 600;2 600;3 600;4 600;
         12 300;13 300;14 300;18 300];
     sysmat=zeros(18,18);
     rhs=zeros(18,1);
    for n=1:size(elements,1)
     [elemat,elevec] =evaluate_stat(nodes(elements(n,:),:),gx2dref(2),gw2dref(2));
     [sysmat,rhs] = assemble(elemat,elevec,sysmat,rhs,elements(n,:));
    end
    [sysmat,rhs] = assignDBC(sysmat,rhs,dbc);
    T=sysmat\rhs;
    r=r+0.01;
end
r=r-0.01
quadplot(nodes,elements,T)
shading interp; 
grid on;
colormap(hot);
colorbar;