McCabe-Thiele | Programa em MATLAB

Prévia do material em texto

Método McCabe-Thiele | Plotagem Gráfica
Programador: Gustavo Meira
Ajuste da Curva de Equilíbrio
Começaremos o métodos desenhando a curva de equilíbrio. Inserimos os pontos coletados de um diagrama 
empírico e ajustamos os dados à essa curva para obter uma equação polinomial aproximada.
% Exemplo: [0,0],[1,3],[2,6],[7,5],[9,4]
% Ajustes usando uma função criada pelo prórpio programador
Coeficientes = polyplotline([0,0],
[0.02,0.14],[0.04,0.13],[0.1,0.42],[0.15,0.52],[0.2,0.58],[0.3,0.67],
[0.4,0.73],[0.45,0.76],[0.5,0.78],[0.55,0.8],[0.6,0.83],[0.65,0.85],
[0.7,0.87],[0.75,0.89],[0.8,0.92],[0.85,0.94],[0.9,0.96],[0.95,0.98],
[1,1],'k');
axis([0 1 0 1]); hold on; plot(0:1:1,0:1:1,'k'); grid minor;
A função usada para plotar o ajuste de curvas é a seguinte:
function p = polyplotline(varargin)
%POLYPLOTLINE plota uma curva ajustada [x,y] inseridos e retorna a função
% Detailed explanation goes here
narginchk(2,inf);
% Inicialização dos valores
1
jj = 0;
especificacao = '';
% Captação dos pontos e das especificações de plotagem
for ii = 1:nargin
 if ischar(varargin{ii})
 especificacao = varargin{ii};
 else
 jj = jj+1;
 x(jj) = varargin{ii}(1);
 y(jj) = varargin{ii}(2);
 end
end
% Plotagem do gráfico
n = input('Qual o grau do polinômio que deseja ajustar?');
p = polyfit(x,y,n);
 if isempty(especificacao)
 x1 = linspace(x(1),x(nargin),100);
 else
 x1 = linspace(x(1),x(nargin-1),100);
 end
f = polyval(p,x1);
 if isempty(especificacao)
 plot(x1,f);
 else
 plot(x1,f,especificacao);
 end
grid on;
Balanço Geral da Coluna
Agora, podemos solucionar a problemática do McCabe-Thiele. O início se dá pelo Balanço Mássico Global:
Após encontrar a variável necessária, podemos seguir para o Balanço Mássico por Componente:
Chamaremos o composto mais volátil de A e faremos o balanço em cima dele, assim como a plotagem do 
gráfico McCabe-Thiele. Sabendo disso, informe a composição de entrada do componente mais volátil:
z_AF = 0.5;
Qual a quantidade do componente mais volátil que você dejesa obter no destilado e no fundo?
% Exemplo: 0.95
y_AD = 0.905;
y_BD = abs(1-y_AD);
x_AB = 0.05;
x_BB = abs(1-x_AB);
Seção de Enriquecimento
2
Podemos começar a plotar as retas de parâmetro do método. Primeiramente, vamos encontrar a razão de 
refluxo r [ ].
r = 1.1653;
Sendo L a vazão de líquido da seção de enriquecimento realizada pelo balanço (V é a vazão de 
vapor). A partir de um balanço de componente, podemos chegar na equação da linha de operação da seç˜ão 
de enriquecimento:
x_enriquecimento = [0 y_AD];
y1 = y_AD/(r+1);
y_enriquecimento = [y1 y_AD];
p_enriquecimento = polyfit(x_enriquecimento,y_enriquecimento,1);
enriquecimento = polyval(p_enriquecimento,linspace(0+0.2,y_AD));
plot([y_AD y_AD],[0 y_AD],'k--');
plot(linspace(0+0.2,y_AD),enriquecimento,'k');
Seção de Alimentação
3
Para definirmos um parâmetro de alimentação, realizamos alguns balanços de energia. A conclusão que 
chegamos é um parâmetro q adimensional que se refere à razão entre a quantidade de calor necessária para 
converter 1 mol de alimentação em vapor saturado e o calor latente molar de vaporização:
q = 0.263;
Lembrando sempre que os valores de entalpia são referentes ao estágio de alimentação! A partir dele, 
Podemos usar a equação a seguir chamada de linha de alimentação para traçar uma reta muito importante:
x_alimentacao = [0 z_AF];
y2 = -z_AF/(q-1);
y_alimentacao = [y2 z_AF];
p_alimetacao = polyfit(x_alimentacao,y_alimentacao,1);
alimetacao = polyval(p_alimetacao,linspace(0,z_AF));
plot([z_AF z_AF],[0 z_AF],'k--');
plot(linspace(0,z_AF),alimetacao,'k');
Seção de Esgotamento
4
A equação que define a linha de operação da seção de esgotamento é a seguinte (sendo o balanço 
):
Podemos, contudo, plotar a linha de esgotamento a partir da fração mássica do componente menos volátil no 
fundo e com a intersecção entre as linhas de operação da seção de enriquecimento e alimentação!
Resolva o sistema abaixo com os valor dos coeficientes respectivos à da alimentação e do 
enriquecimento para obter o ponto p.
disp(p_alimetacao);
 -0.3569 0.6784
disp(p_enriquecimento);
 0.5382 0.4180
p1 = 0.2909;
x_esgotamento = [x_AB p1];
p2 = polyval(p_enriquecimento,p1);
y_esgotamento = [x_AB p2];
p_esgotamento = polyfit(x_esgotamento,y_esgotamento,1);
esgotamento = polyval(p_esgotamento,linspace(x_AB,p1+0.1));
plot([x_AB x_AB],[0 x_AB],'k--');
plot(linspace(x_AB,p1+0.1),esgotamento,'k');
5
Número Ideal de Pratos
O número ideal de pratos é identificado a partir das iterações feitas a partir das retas horizontais e verticais 
traçadas, respectivamente e sucessivamente, a partir do ponto . As iterações só param quando o ponto 
é alcançado e ultrapassado.
moving_x = y_AD;
pratos = 0;
while moving_x > x_AB
 % Plotagem da reta horizontal
 if moving_x >= p1
 y1 = polyval(p_enriquecimento,moving_x);
 x1 = @(x)polyval(p_enriquecimento,x)-y1;
 x1 = fzero(x1,moving_x);
 x2 = @(x)polyval(Coeficientes,x)-y1;
 x2 = fzero(x2,moving_x);
 x_plot = [x2 x1];
 y_plot = [y1 y1];
 p_horizontal = polyfit(x_plot,y_plot,1);
 horizontal = polyval(p_horizontal,linspace(x2,x1));
 plot(linspace(x2,x1),horizontal,'Color',[0.9290 0.6940 0.1250]);
 else
 y1 = polyval(p_esgotamento,moving_x);
 x1 = @(x)polyval(p_esgotamento,x)-y1;
 x1 = fzero(x1,moving_x);
6
 x2 = @(x)polyval(Coeficientes,x)-y1;
 x2 = fzero(x2,moving_x);
 x_plot = [x2 x1];
 y_plot = [y1 y1];
 p_horizontal = polyfit(x_plot,y_plot,1);
 horizontal = polyval(p_horizontal,linspace(x2,x1));
 plot(linspace(x2,x1),horizontal,'Color',[0.9290 0.6940 0.1250]);
 end
 moving_x = x2;
 pratos = pratos +1;
 % Plotagem da reta vertical
 if moving_x >= p1
 x_plot = [x2 x2];
 y2 = polyval(p_enriquecimento,x2);
 y_plot = [y1 y2];
 plot(x_plot,y_plot,'Color',[0.9290 0.6940 0.1250]);
 else
 x_plot = [x2 x2];
 y2 = polyval(p_esgotamento,x2);
 y_plot = [y1 y2];
 plot(x_plot,y_plot,'Color',[0.9290 0.6940 0.1250]);
 end
end
title('McCabe-Thiele'); hold off
xlabel('x_{líquido}'); ylabel('y_{vapor}');
legend('','','','','','','','',['Nº Pratos = ', 
num2str(pratos)],'Location','northwest');
7
8