Tesina Modelli e metodi numerici trasporto refrigerato

Formattazione del testo con tag HTML

'T [°C]')title(color,% [Tx, Ty] = Me.gradient(T_sol);% figure% Me.quiver(Tx, Ty)% Soluzione del transitorio% Datitfin = 60*60*5; % Intervallo di tempo in analisidt = 100; % Passo temporalent = tfin/dt; % Numero di passi temporali% Tempo di variazione della T dell'aria neltmax = 600;furgone % Vettore della temperatura nei dofT_dof = T_sol(3:end) + 273.15;T_puntuale = zeros(nt, 4); % Inizializzazione della matrice in cuiverrà salvato l'andamento della T in 4 punti diversiT_aria = zeros(nt, 1); % Inizializzazione del vettore in cui verràsalvato l'andamento della T dell'aria nel furgonestabEE = 2/abs(eigs(D, M, 1)); % Condizione di stabilità di EE% Matrice costante per il metodo di EIE = (M + dt*D);for k = 1:ntT_aria(k) = T_int - 273.15; % Controllo T aria nel furgoneif dt*k <= tmax % Se la T dell'aria nel furgone cambia,cambiano anche le condizioni di Robin, quindi cambia il vettore bT_int = T_int + dt*k*(-6 + 273.15 -

T_int)/(tmax); % Aggiornamento Taria nel furgone
% Aggiornamento delle condizioni di Robin
Pallet1.Borders.Bc([1, 5, 7]) = boundaries.robin(alfa_18_or1,(alfa_18_or1)*T_int);
Pallet1.Borders.Bc([2, 4, 6, 8]) = boundaries.robin(alfa_18_vert1,(alfa_18_vert1)*T_int);
Pallet2.Borders.Bc([2, 4]) = boundaries.robin(alfa_18_vert2,(alfa_18_vert2)*T_int);
Pallet2.Borders.Bc(3) = boundaries.robin(alfa_18_or2,(alfa_18_or2)*T_int);
Canale.Borders.Bc([3, 5, 7]) = boundaries.robin(alfa_28,(alfa_28)*T_int);
Furgone = [Pallet1, Pallet2, Canale]; % Aggiornamento di Me
Me = mesh2D(Furgone, 0.01); % Aggiornamento di b
[D,b] = BuildMatTruckStab(Me,f);
end % Calcolo della soluzione al
T_dof = E \(M*T_dof + dt*b); % passo temporale dt*k con EI
% Costruzione del vettore
T = [-28; -28; T_dof - 273.15]; % soluzione in tutti i nodi
Me.draw(T);
shading interp
title(['t=',num2str(k*dt)])
pause(0.005)
hold off
% Registrazione della T nell'istante dt*k in quattro punti notevoli
% Punta dx Pallet 1
T_puntuale(k, 1) = T(8);
T_puntuale(k,

Transport-Diffusion-Reaction problem with non homogeneous Dirichlet B.C.s

Input:

• Me: a Mesh2D object
• f: MATLAB function of (x,y) which returns the values of the external source. Default: constant value=4

Output:

• D: diffusion+transport+reaction matrix
• b: constant terms vector

Check inputs

if nargin < 2
    f = @(x,y) 4*ones(size(x));
end

For clarity, call some properties of Me with shorter names

V = Me.Triangles.Vertices;
Areas = Me.Triangles.Areas;
CenterOfMass = Me.Triangles.CenterOfMass;
Nodes = Me.Nodes;
Dof = Me.Nodes.Dof;

Number of internal nodes: we know that the N unknown nodes are numbered from 1 to N in Me.UnknownNodes; the maximum is therefore the number of unknown (degrees of freedom)

numDof = max(Dof);

Vectors preallocation: instead of allocating the (sparse) diffusion matrix, we save the rows, columns, and values corresponding to each contribution; at the end, we'll call sparse(...) to obtain the diffusion matrix

b = zeros(numDof,1);
row = zeros(Me.MatrixContributions,1);
col = ...

zeros(Me.MatrixContributions,1);
d=zeros(Me.MatrixContributions,1);
pos=1; 
%we start from the element in position 1, we'll increase this index
%everytime we add an entry
%we evaluate the external force in the center of mass of this triangle
force = f(CenterOfMass.X,CenterOfMass.Y);
%evaluate the value of the coefficient in front of the Laplace operator
mu=Me.mu;
%evaluate the value of the coefficient in front of the divergence operator
%NB. beta is a two-coefficient vector, indicating the speed in the x
%and y directions
beta=Me.beta;
rho=Me.rho;
%main loop on each triangle
for e=1:size(V,1)
    Dx(1) = Nodes.X(V(e,3)) - Nodes.X(V(e,2));
    Dx(2) = Nodes.X(V(e,1)) - Nodes.X(V(e,3));
    Dx(3) = Nodes.X(V(e,2)) - Nodes.X(V(e,1));
    Dy(1) = Nodes.Y(V(e,3)) - Nodes.Y(V(e,2));
    Dy(2) = Nodes.Y(V(e,1)) - Nodes.Y(V(e,3));
    Dy(3) = Nodes.Y(V(e,2)) - Nodes.Y(V(e,1));
    %for each vertex of this triangle
    for ni=1:3
        %look at the "unknown" numbering: if the node is positive, it
        %corresponds to a degree of freedom of the

problemii = Dof(V(e,ni));%is it unknown?
if ii > 0%yes it is! second loop on the vertices
    for nj=1:3
        dtmp=mu(e)*(Dy(ni)*Dy(nj)+Dx(ni)*Dx(nj))/(4.0*Areas(e)) + ...(beta(e,1)*Dy(nj)-beta(e,2)*Dx(nj))*rho(e)/6;
        jj = Dof(V(e,nj));%%is it unknown as well?
        if jj > 0
            row(pos)=ii;
            col(pos)=jj;
            d(pos)=dtmp;
            pos=pos+1;
        else
            val=Me.BC.DirichletNodes(-jj,2);
            b(ii)=b(ii)-dtmp*val;
        end
    end
    %build the constant terms vector adding the external%contribution
    b(ii) = b(ii) + Areas(e)*force(e)/3.0;
end
end
end
%% Robin B.C.s
Edges=Me.Edges;
Robin=Me.BC.RobinEdges;
for k=1:size(Robin,1)
    Node1=Edges(Robin(k,1),1);
    Node2=Edges(Robin(k,1),2);
    dx=Nodes.X(Node1)-Nodes.X(Node2);
    dy=Nodes.Y(Node1)-Nodes.Y(Node2);
    dist=sqrt(dx*dx+dy*dy);
    ii1=Dof(Node1);
    ii2=Dof(Node2);
    g=Robin(k,3);
    h=Robin(k,2);
    if ii1>0 && ii2<0 %ii1 is unknown, ii2 is known
        Dirichlet=Me.BC.DirichletNodes(-ii2,2);
        b(ii1)=b(ii1)+g/2*dist-h*dist/6*Dirichlet;
        row(pos)=ii1;
        col(pos)=ii1;
        d(pos)=h*dist/3;
        pos=pos+1;
        %D(ii1,ii1)=D(ii1,ii1)+h*dist/3;
    elseif ii1<0

<code>&amp;&amp; ii2>0 %ii1 is known, ii2 is unknown
Dirichlet=Me.BC.DirichletNodes(-ii1,2);
b(ii2)=b(ii2)+g/2*dist-h*dist/6*Dirichlet;
row(pos)=ii2;
col(pos)=ii2;
d(pos)=h*dist/3;
pos=pos+1;
%D(ii2,ii2)=D(ii2,ii2)+h*dist/3;
else %both are unknwon
b(ii1)=b(ii1)+g/2*dist;
b(ii2)=b(ii2)+g/2*dist;
row(pos:pos+3)=[ii1;ii2;ii1;ii2];
col(pos:pos+3)=[ii1;ii2;ii2;ii1];
d(pos:pos+3)=[2;2;1;1]*h*dist/6;
pos=pos+4;
end
end
%% Neumann B.C.s
Edges=Me.Edges;
Neumann=Me.BC.NeumannEdges;
for k=1:size(Neumann,1)
Node1=Edges(Neumann(k),1);
Node2=Edges(Neumann(k),2);
dx=Nodes.X(Node1)-Nodes.X(Node2);
dy=Nodes.Y(Node1)-Nodes.Y(Node2);
dist=sqrt(dx*dx+dy*dy);
g=Me.BC.NeumannEdges(k,2);
if Dof(Node1)>0
b(Dof(Node1))=b(Dof(Node1))+g/2*dist;
end
if Dof(Node2)>0
b(Dof(Node2))=b(Dof(Node2))+g/2*dist;
end
end
end
%% assemble the stiffness matrix
DD = sparse(row(1:pos-1), col(1:pos-1), d(1:pos-1), numDof, numDof);
end FUNZIONI BUILD MAT TRUCK STAB:
function [D,b] = BuildMatTruckStab(Me,f)
%Assemble the matrix D and the vector b of the</code>

Transport-Diffusion-Reaction problem with non homogeneous Dirichlet B.C.s with artificial diffusion

Input:

• Me: a Mesh2D object
• f: MATLAB function of (x,y) which returns the values of the external source. Default: constant value=4

Output:

• D: diffusion+transport+reaction matrix
• b: constant terms vector

Check inputs

if nargin < 2
    f = @(x,y) 4*ones(size(x));
end

For clarity, call some properties of Me with shorter names

V = Me.Triangles.Vertices;
Areas = Me.Triangles.Areas;
CenterOfMass = Me.Triangles.CenterOfMass;
Nodes = Me.Nodes;
Dof = Me.Nodes.Dof;

Number of internal nodes: we know that the N unknown nodes are numbered from 1 to N in Me.UnknownNodes; the maximum is therefore the number of unknown (degrees of freedom)

numDof = max(Dof);

Vectors preallocation: instead of allocating the (sparse) diffusion matrix, we save the rows, columns and values corresponding to each contribution; at the end, we'll call sparse(...) to obtain the diffusion matrix

b = zeros(numDof,1);
row = 

zeros(Me.MatrixContributions,1);

node is positive, it%corresponds to a degree of freedom of the problemii = Dof(V(e,ni));%is it unknown?if ii > 0%yes it is! second loop on the verticesfor nj=1:3if Peclet(e)>1muStab=mu(e)*(1+Peclet(e));else muStab=mu(e);enddtmp=muStab*(Dy(ni)*Dy(nj)+Dx(ni)*Dx(nj))/(4.0*Areas(e)) + ...(beta(e,1)*Dy(nj)-beta(e,2)*Dx(nj))*rho(e)/6;jj = Dof(V(e,nj));%%is it unknown as well?if jj > 0row(pos)=ii;col(pos)=jj;d(pos)=dtmp;pos=pos+1;else val=Me.BC.DirichletNodes(-jj,2);b(ii)=b(ii)-dtmp*val;endend%build the constant terms vector adding the external%contributionb(ii) = b(ii) + Areas(e)*force(e)/3.0;endendend%% Robin B.C.sEdges=Me.Edges;Robin=Me.BC.RobinEdges;for k=1:size(Robin,1)Node1=Edges(Robin(k,1),1);Node2=Edges(Robin(k,1),2);dx=Nodes.X(Node1)-Nodes.X(Node2);dy=Nodes.Y(Node1)-Nodes.Y(Node2);dist=sqrt(dx*dx+dy*dy);ii1=Dof(Node1);ii2=Dof(Node2);g=Robin(k,3);