thermal_stat.m

study of rod with finite element method

% Study of the thermal flow between a concrete wall and a concrete slab
%
clear all;   clc;   close all;
%
% Geometrical data and generation of the mesh
%
    nodes = [];   elements = [];   n_elem = 0;
% First super-element
    i1 = 1;	% number of the node number 1 in local co-ordinates
    subd_x = [13 9 7 6 5];   subd_y  = [5 5];
    n1 = length(subd_x);   n2 = length(subd_y);   % number of subdivisions in the x and y directions
    xl = -65;
    for i = 1:n1
        xl = [xl xl(2*i-1)+subd_x(i)/2*[1 2]];
    end
    yl = 0;
    for i = 1:n2
        yl = [yl yl(2*i-1)+subd_y(i)/2*[1 2]];
    end
    N1 = 2*n1+1;   N2 = 2*n2+1;   n_elem = n_elem+n1*n2;
    nodes = [xl(1)*ones(N2,1) yl'];
    for i = 1:n1
        nodes = [ nodes ; [xl(2*i)*ones(n2+1,1) yl(1:2:N2)'] ; [xl(2*i+1)*ones(N2,1) yl'] ];
    end
    elements = [ elements ; i1+[0 N2 3*n2+1+[1:3] N2+1 [2 1]] ];
    [i2,i3] = size(elements);
    for i = i2+[1:n2-1]
        elements(i,:) = elements(i-1,:)+2;
        elements(i,2) = elements(i-1,2)+1;
        elements(i,6) = elements(i-1,6)+1;
    end
    for i = 1:n1-1
        elements = [elements ; (3*n2+2)*i+elements(i2+[0:n2-1],:)];
    end
% Second super-element
    i1 = length(nodes)+1-N2;	% number of the node number 1 in local co-ordinates
    subd_x = [5 5 5 5 5 5 5 5 5 5];   subd_y  = [5 5 5 5 8 10 12];
    n1 = length(subd_x);   n2 = length(subd_y);   % number of subdivisions in the x and y directions
    xl = nodes(i1,1);
    for i = 1:n1
        xl = [xl xl(2*i-1)+subd_x(i)/2*[1 2]];
    end
    yl = nodes(i1,2);
    for i = 1:n2
        yl = [yl yl(2*i-1)+subd_y(i)/2*[1 2]];
    end
    N1 = 2*n1+1;   N2 = 2*n2+1;   n_elem = n_elem+n1*n2;
    nodes(i1:i1+N2-1,:) = [xl(1)*ones(N2,1) yl'];
    for i = 1:n1
        nodes = [ nodes ; [xl(2*i)*ones(n2+1,1) yl(1:2:N2)'] ; [xl(2*i+1)*ones(N2,1) yl'] ];
    end
    elements = [ elements ; i1+[0 N2 3*n2+1+[1:3] N2+1 [2 1]] ];
    [i2,i3] = size(elements);
    for i = i2+[1:n2-1]
        elements(i,:) = elements(i-1,:)+2;
        elements(i,2) = elements(i-1,2)+1;
        elements(i,6) = elements(i-1,6)+1;
    end
    for i = 1:n1-1
        elements = [elements ; (3*n2+2)*i+elements(i2+[0:n2-1],:)];
    end
% Third super-element
    i1 = length(nodes)+1-N2;	% number of the node number 1 in local co-ordinates
    i2 = i1+10;   % fictive number of the node number 1 in local co-ordinates
    ielem = n_elem;
    subd_x = [5 6 7 9 13];   subd_y  = [5 5];
    n1 = length(subd_x);   n2 = length(subd_y);   % number of subdivisions in the x and y directions
    xl = nodes(i1,1);
    for i = 1:n1
        xl = [xl xl(2*i-1)+subd_x(i)/2*[1 2]];
    end
    yl = nodes(i1,2);
    for i = 1:n2
        yl = [yl yl(2*i-1)+subd_y(i)/2*[1 2]];
    end
    N1 = 2*n1+1;   N2 = 2*n2+1;   n_elem = n_elem+n1*n2;
    for i = 1:n1
        nodes = [ nodes ; [xl(2*i)*ones(n2+1,1) yl(1:2:N2)'] ; [xl(2*i+1)*ones(N2,1) yl'] ];
    end
    elements = [ elements ; i2+[0 N2 3*n2+1+[1:3] N2+1 [2 1]] ];
    [i2,i3] = size(elements);
    for i = i2+[1:n2-1]
        elements(i,:) = elements(i-1,:)+2;
        elements(i,2) = elements(i-1,2)+1;
        elements(i,6) = elements(i-1,6)+1;
    end
    for i = 1:n1-1
        elements = [elements ; (3*n2+2)*i+elements(i2+[0:n2-1],:)];
    end
    for i = 1:n2
        elements(ielem+i,[1 8 7]) = i1+2*(i-1)+[0:2];
    end
    nnodes = length(nodes);
%
% Plot of the nodes and elements
%
    figure;   hold on;
    plot(nodes(:,1),nodes(:,2),'ro','LineWidth',2);
    xlabel('abscis x [cm]');
    ylabel('distance to symmetry axis y [cm]');
    for i = 1:n_elem
        n = [elements(i,:) elements(i,1)];
        plot(nodes(n,1),nodes(n,2),'k-','LineWidth',2);
    end
    legend('nodes','elements',0);
    grid on;   axis image;
%     disp(' ');   disp(' ');   disp(' ');
%     answer = input(['          Do you agree with the plotted mesh ?',...
%         '\n','                            If you agree, pull    y   :           '],'s');
%     if not(answer == 'y' | answer == 'Y')   break;   end
%
% Imposed nodal temperatures
%
    thetap(nnodes) = 0;
    nt1 = [];   nt2 = [];   nt3 = [];
    for i = 1: nnodes
        if abs(nodes(i,1)) < .0001
            nt1 = [nt1 i];
        end
        if nodes(i,1) < -24.9999 & nodes(i,2) > 9.9999
            nt2 = [nt2 i];
        end
        if nodes(i,1) > 24.9999 & nodes(i,2) > 9.9999
            nt3 = [nt3 i];
        end
    end
    A = [nt3' -nodes(nt3',2)];   B = sortrows(A,2);   nt3 = B(:,1);
    thetap(nt1) = -5*ones(size(nt1));
    thetap(nt2) = 15*ones(size(nt2));
    thetap(nt3) = 15*ones(size(nt2));
    figure;   hold on;
    plot(nodes(:,1),nodes(:,2),'ro','LineWidth',2);
    xlabel('abscis x [m]');
    ylabel('distance to symmetry axis y [m]');
    for i = 1:n_elem
        n = [elements(i,:) elements(i,1)];
        plot(nodes(n,1),nodes(n,2),'k-','LineWidth',2);
    end
    grid on;   axis image;
    plot(nodes(nt1,1),nodes(nt1,2),'b-',nodes(nt2,1),nodes(nt2,2),'r-',nodes(nt3,1),...
        nodes(nt3,2),'r-','LineWidth',4);
    text(1,25,'T = -5?C','Color','blue');
    text(27,13,'T = 15?C','Color','red');
    text(-33,13,'T = 15?C','Color','red');
    hold off;
%
% Inserting the values of the thermal conductivity coefficient
%
    lambda = ones(1,n_elem);   lambda([27:31 60:66]) = .05;
%
% Calculation of the form function
%
xe = [-1 0 1 1 1 0 -1 -1]';
ye = [-1 -1 -1 0 1 1 1 0]';
    aij = [ones(8,1) xe ye xe.*ye xe.^2 ye.^2 (xe.^2).*ye xe.*(ye.^2)];
    cij = inv(aij);
%
% Numerical calculation of the element stiffness matrix Ke (first case)
%
    xi = sqrt(.6)*[-1 -1 -1 0 0 0 1 1 1];
    yi = sqrt(.6)*[-1 0 1 -1 0 1 -1 0 1];
    wi = [25 40 25 40 64 40 25 40 25]/81;
    K(nnodes,nnodes) = 0;
    for i = 1:n_elem
        n = elements(i,:);   xxe = nodes(n,:);
        Ke = [];   Ke(8,8) = 0;
        for ip = 1:length(xi)
            D = [0 1 0 yi(ip) 2*xi(ip) 0 2*xi(ip)*yi(ip) yi(ip)^2 ;
                0 0 1 xi(ip) 0 2*yi(ip) xi(ip)^2 2*xi(ip)*yi(ip) ]*cij;
            J = D*xxe;   detJ = det(J);   grad = inv(J)*D;
            Ke = Ke + wi(ip)*det(J)*lambda(i)*grad'*grad;
        end
        K(n,n) = K(n,n)+Ke;
    end
    Kmod = K;
    for i = 1:nnodes
        if thetap(i)~=0
            Kmod(i,:) =  [zeros(1,i-1) 1 zeros(1,nnodes-i)];
        end
    end
    theta = Kmod\thetap';
    r1 = [1:170]';
    r2 = [];
    start = nnodes;
    for i = 1:5
        r2 = [r2 ; start-8*(i-1)+[-4:0 -7:-5]' ];
    end
    start = nnodes-length(r2);
    for i = 1:5
        r2 = [r2 ; start-23*(i-1)+[-14:0 -22:-15]' ];
    end
    r2 = [r2 ; [156:170]' ];
    nodal_temperature = [r1 theta(r1) theta(r2)];
    disp(' ');   disp(' ');   disp(' ');
    disp('Node Nr     Not isolated      Well isolated')
    for i = 1:170
        disp(sprintf('%10i %20.4f %20.4f',nodal_temperature(i,:)));
    end
    disp(' ');   disp(' ');   disp(' ');
%
% Plot of the band structure
%
    figure;   spy([Kmod zeros(nnodes,2) thetap']);
%
% Calculation of interpolation value for the different plots
%
    X = -1:.2:1;
    Y = -1:.2:1;
    [x,y] = meshgrid(X,Y);
    i = 1;
    ne1 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 2;
    ne2 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 3;
    ne3 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 4;
    ne4 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 5;
    ne5 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 6;
    ne6 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 7;
    ne7 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
    i = 8;
    ne8 = ones(size(x))*cij(1,i)+x*cij(2,i)+y*cij(3,i)+x.*y*cij(4,i)+x.^2*cij(5,i)+...
        y.^2*cij(6,i)+x.^2.*y*cij(7,i)+x.*y.^2*cij(8,i);
%
% Plot of the temperature distribution in the first case
%
    figure;   hold on;
    xnodes = nodes(:,1);   ynodes = nodes(:,2);
    for i = 1:n_elem
        n = elements(i,:);
        xp = ne1*xnodes(n(1))+ne2*xnodes(n(2))+ne3*xnodes(n(3))+ne4*xnodes(n(4))+...
            ne5*xnodes(n(5))+ne6*xnodes(n(6))+ne7*xnodes(n(7))+ne8*xnodes(n(8));
        yp = ne1*ynodes(n(1))+ne2*ynodes(n(2))+ne3*ynodes(n(3))+ne4*ynodes(n(4))+...
            ne5*ynodes(n(5))+ne6*ynodes(n(6))+ne7*ynodes(n(7))+ne8*ynodes(n(8));
        tp = ne1*theta(n(1))+ne2*theta(n(2))+ne3*theta(n(3))+ne4*theta(n(4))+...
            ne5*theta(n(5))+ne6*theta(n(6))+ne7*theta(n(7))+ne8*theta(n(8));
        plot3(nodes(n,1),nodes(n,2),16*ones(size(n)),'k-','LineWidth',1);
        [c1,h1] = contour('v6',xp,yp,tp,-5:.5:15);   set(h1,'LineWidth',2);
    end
title('Wall with thermal bridge        7 8-noded quadrangular elements           Perfectly isolated wall');
    plot3(xnodes,ynodes,16*ones(size(xnodes)),'ko',...
        xnodes,ynodes,16*ones(size(xnodes)),'k.','LineWidth',2);
    axis image;   colorbar('hor');   hold off;
%
% Plot of the temperature variation
%
% Along the horizontal symmetry axis
    figure;   hold on;   axis([-66 1 -6 16]);
    nx1 = [1 6 9 14 17 22 25 30 33 38 41 56 64 79 87 102 110 125 133 148 156];
    x1 = xnodes(nx1);   t11 = nodal_temperature(nx1,2);   t12 = nodal_temperature(nx1,3);
    h = plot(x1,t11,'r-',x1,t12,'b-','LineWidth',2);
    hold on;
% Along the top of the slab
    nx2 = [5 8 13 16 21 24 29 32 37 40 45 58 68 81 91 104 114 127 137 150 160];
    x2 = xnodes(nx2);   t21 = nodal_temperature(nx2,2);   t22 = nodal_temperature(nx2,3);
    plot(x2,t21,'m:',x1,t22,'c:',-[0 10 15 25 65],[-15 -10 40 45 45]/3,'k*','LineWidth',2);
    grid on;
    title('Temperature variation - 7 8-noded quadrangular elements');
    xlabel('abscis  [cm]');   ylabel('temperature  [?C]');
    legend('Symmetry line - not isolated','Symmetry line - isolated',...
        'Top of the slab - not isolated','Top of the slab - isolated',...
        'Theoretical values - isolated',0);
    axis tight;   hold off;

% Calculation of the total thermal losses
%
    xq = xi(7:9);   yq = yi(7:9);   wq = [5 8 5]/9;
    qx = [0 0];
    for i = 39:45
       n = elements(i,:);
       for iq = 1:length(xq)
           D = [0 1 0 yq(iq) 2*xq(iq) 0 2*xq(iq)*yq(iq) yq(iq)^2 ;
               0 0 1 xq(iq) 0 2*yq(iq) xq(iq)^2 2*xq(iq)*yq(iq) ]*cij;
           J = D*nodes(n,:);   detJ = J(2,2);
           grad = inv(J)*D*[lambda(i)*theta(n,:)];
           qx(1) = qx(1)+wq(iq)*detJ*grad(1,:);
       end
    end
    xq = xi(1:3);   yq = yi(1:3);   wq = [5 8 5]/9;
    for i = 46:52
        n = elements(i,:);
        for iq = 1:length(xq)
            D = [0 1 0 yq(iq) 2*xq(iq) 0 2*xq(iq)*yq(iq) yq(iq)^2 ;
                0 0 1 xq(iq) 0 2*yq(iq) xq(iq)^2 2*xq(iq)*yq(iq) ]*cij;
           J = D*nodes(n,:);   detJ = J(2,2);
           grad = inv(J)*D*[lambda(i)*theta(n,:)];
            qx(2) = qx(2)+wq(iq)*detJ*grad(1,:);
        end
    end
    heat_flow = qx