bottennoder för referens är 1 6 12 19 27 35 43 51 59 67 75 83 91 98 104

clear all

Coord = [ 0 0; 0 5; 0 10; 0 15; 0 20; 20/3 0; 20/3 5; 20/3 10; 20/3 15; 20/3 20; 20/3 25; 40/3 0; 40/3 5; 40/3 10; 40/3 15; 40/3 20; 40/3 25; 40/3 30; 20 0; 20 5; 20 10; 20 15; 20 20; 20 25; 20 30; 20 35; 25 0; 25 5; 25 10; 25 15; 25 20; 25 25; 25 30; 25 35; 30 0; 30 5; 30 10; 30 15; 30 20; 30 25; 30 30; 30 35; 35 0; 35 5; 35 10; 35 15; 35 20; 35 25; 35 30; 35 35; 40 0; 40 5; 40 10; 40 15; 40 20; 40 25; 40 30; 40 35; 45 0; 45 5; 45 10; 45 15; 45 20; 45 25; 45 30; 45 35; 50 0; 50 5; 50 10; 50 15; 50 20; 50 25; 50 30; 50 35; 55 0; 55 5; 55 10; 55 15; 55 20; 55 25; 55 30; 55 35; 60 0; 60 5; 60 10; 60 15; 60 20; 60 25; 60 30; 60 35; 200/3 0; 200/3 5; 200/3 10; 200/3 15; 200/3 20; 200/3 25; 200/3 30; 220/3 0; 220/3 5; 220/3 10; 220/3 15; 220/3 20; 220/3 25; 240/3 0; 240/3 5; 240/3 10; 240/3 15; 240/3 20 ];

plot(Coord(:,1), Coord(:,2),'x')

%konstanter
kasf = 30;
ksteel = 45;
casf = 1200;
eq = 7*10^5;
Tg = 25; %Celcius
Tinf = 25; %Celcius

ep = 30; %mm
D1 = kasf*eye(2);
D2 = ksteel*eye(2);

c1 = 1/6*[ 2 0 1; 0 0 0; 1 0 2 ];
c2 = 1/6*[ 2 1 0; 1 2 0; 0 0 0 ];

L1 = 5;
L2 = sqrt(5^2 + (20/3)^2);

KE = zeros(3,3,178);
Kce = zeros(3,3,178);
fe = zeros(3,1);
Edof = zeros(4,178);
Dof = (1:108)';
K = zeros(108,108);
%Kc = zeros(108,108);
f = zeros(108, 1);
f(1) = eq;
f(6) = eq;
f(12) = eq;
f(19) = eq;
f(27) = eq;
f(35) = eq;
f(43) = eq;
f(51) = eq;
f(59) = eq;
f(67) = eq;
f(75) = eq;
f(83) = eq;
f(91) = eq;
f(98) = eq;
f(104) = eq;

l = 0;
for j=1:49,
    i = Coord(j,1);

    if (i == 0)
        nodejump = 5;
    elseif (i == 20/3)
        nodejump = 6;
    elseif (i == 40/3)
        nodejump = 7;
    else
        nodejump = 8;
    end

            if (Coord(j,2) < Coord(j+nodejump+1,2))
                 l = l + 1;

                x1 = Coord(j,1);
                y1 = Coord(j,2);

                x2 = Coord(j+nodejump,1);
                y2 = Coord(j+nodejump,2);

                x3 = Coord(j+nodejump+1,1);
                y3 = Coord(j+nodejump+1,2);

                Edof(:,l) = [l j j+nodejump j+nodejump+1 ];

                ex = [x1 x2 x3];
                ey = [y1 y2 y3];

                if ((j == 38) || (j == 39) || (j == 46) || (j == 47))
                    KE(:,:,l) = flw2te(ex,ey,ep,D2);
                else
                    KE(:,:,l) = flw2te(ex,ey,ep,D1);
                end
            end

            if (Coord(j,2) < Coord(j+1,2))
                l = l + 1;

                x1 = Coord(j,1);
                y1 = Coord(j,2);

                x2 = Coord(j+nodejump+1,1);
                y2 = Coord(j+nodejump+1,2);

                x3 = Coord(j+1,1);
                y3 = Coord(j+1,2);



                Edof(:,l) = [l j j+nodejump+1 j+1 ];

                ex = [x1 x2 x3];
                ey = [y1 y2 y3];

                if ((j == 38) || (j == 39) || (j == 46) || (j == 47))
                    KE(:,:,l) = flw2te(ex,ey,ep,D2);
                else
                    KE(:,:,l) = flw2te(ex,ey,ep,D1);
                end
            end
end

l = 89;
for j=51:102,
    i = Coord(j,1);

    if (i == 200/3)
        nodejump = 7;
    elseif (i == 220/3)
        nodejump = 6;
    else
        nodejump = 8;
    end

            if (Coord(j,2) < Coord(j+1,2))
                 l = l + 1;

                x1 = Coord(j+nodejump,1);
                y1 = Coord(j+nodejump,2);

                x2 = Coord(j+1,1);
                y2 = Coord(j+1,2);

                x3 = Coord(j,1);
                y3 = Coord(j,2);



                Edof(:,l) = [l j+nodejump j+1 j ];

                ex = [x1 x2 x3];
                ey = [y1 y2 y3];

                if ((j == 54) || (j == 55) || (j == 62) || (j == 63))
                   KE(:,:,l) = flw2te(ex,ey,ep,D2);
                else
                   KE(:,:,l) = flw2te(ex,ey,ep,D1);
                end
            end

            if ((j < 102) && Coord(j,2) < Coord(j+nodejump+1,2))
                l = l + 1;

                x1 = Coord(j+nodejump,1);
                y1 = Coord(j+nodejump,2);

                x2 = Coord(j+nodejump+1,1);
                y2 = Coord(j+nodejump+1,2);

                x3 = Coord(j+1,1);
                y3 = Coord(j+1,2);




                Edof(:,l) = [l j+nodejump j+nodejump+1 j+1 ];


                ex = [x1 x2 x3];
                ey = [y1 y2 y3];


                if ((j == 54) || (j == 55) || (j == 62) || (j == 63))
                   KE(:,:,l) = flw2te(ex,ey,ep,D2);
                else
                   KE(:,:,l) = flw2te(ex,ey,ep,D1);
                end
            end
end

Edof = Edof';
[Ex,Ey] = coordxtr(Edof,Coord,Dof,3);
eldraw2(Ex,Ey,[1,2,2])

for j=1:178,
    K=assem(Edof,K,KE(:,:,j));
end

%Boundary cond.
% Tg = 300;
% bc = [(1:108)' zeros(108,1)];
% bc(:,2) = Tinf;
% bc(1,2) = Tg;
% bc(6,2) = Tg;
% bc(12,2) = Tg;
% bc(19,2) = Tg;
% bc(27,2) = Tg;
% bc(35,2) = Tg;
% bc(43,2) = Tg;
% bc(51,2) = Tg;
% bc(59,2) = Tg;
% bc(67,2) = Tg;
% bc(75,2) = Tg;
% bc(83,2) = Tg;
% bc(91,2) = Tg;
% bc(98,2) = Tg;
% bc(104,2) = Tg;

bc = [
%         2 Tinf;
%         3 Tinf;
%         4 Tinf;
%         5 Tinf;
%         11 Tinf;
%         18 Tinf;

        26 Tg;
        34 Tg;
        42 Tg;
        50 Tg;
        58 Tg;
        66 Tg;
        74 Tg;
        82 Tg;
        90 Tg ];

%         97 Tinf;
%         103 Tinf;
%         105 Tinf;
%         106 Tinf;
%         107 Tinf;
%         108 Tinf];

% dt = 0.5;
% T = 60
% ip = [dt T alfa
[a,r]=solveq(K,f,bc);

Ed=extract(Edof,a);
for j=1:178
    if ((j == 38) || (j == 39) || (j == 46) || (j == 47) || (j == 54) || (j == 55) || (j == 62) || (j == 63))
        Es(j,:)=flw2ts(Ex(j,:),Ey(j,:),D2,Ed(j,:));
    else
        Es(j,:)=flw2ts(Ex(j,:),Ey(j,:),D1,Ed(j,:));
    end
end

sfac=scalfact2(Ex,Ey,Es,0.5);
eldraw2(Ex,Ey,[1,3,0]);
elflux2(Ex,Ey,Es,[1,4],sfac);
pltscalb2(sfac,[2e-2 0.06 0.01],4);
pause; clf;
eldraw2(Ex,Ey,[1,3,0]);
eliso2(Ex,Ey,Ed,5,[1,4]);

colormap('jet')
fill(Ex',Ey',Ed')
axis equal