Untitled

function main
    
    clear all;
    close all;
    clc;

    % sistem parametreleri
    global F0 V CA0 T0 Cp rho k0 E R dH Th0 Cph rhoh Fh0
    F0 = 1;             % Besleme Akış Hızı [m³/dak]
    V = 5;              % Reaktör Hacmi [m³]
    CA0 = 5;           % Besleme Konsantrasyonu [kmol/m³]
    T0 = 300;          % Besleme Sıcaklığı [K]
    Cp = 1;            % Isı Kapasitesi [kcal/kg.K]
    rho = 1e3;         % Yoğunluk [kg/m³]
    k0 = 1e7;          % Reaksiyonun Hız Sabiti [dak⁻¹]
    E = 1.104e4;       % Aktivasyon Enerjisi [kcal/kmol]
    R = 1.987;         % Gaz Sabiti [kcal/kmol.K]
    dH = 1.1e5;        % Reaksiyon ısısı [kcal/kmol]
    Th0 = 363;         % Ceket Giriş Sıcaklığı [K]
    Cph = 1;           % Ceket Sıvısının Isı Kapasitesi [kcal/kg.K]
    rhoh = 1e3;        % Ceket Sıvısı Yoğunluğu [kg/m³]
    Fh0 = 20;          % Başlangıç Ceket Sıvı Akış Hızı [m³/dak]

    % Zaman aralığı ve başlangıç koşulları
    t_initial = linspace(0, 20, 1000);  
    t_final = linspace(20, 40, 1000);   
    y0 = [CA0; T0];                    

    % Başlangıç akış hızı için ODEs çöz
    [t1, y1] = ode15s(@CSTR_ODEs, t_initial, y0);
    
    % Akış hızını değiştirip tekrar çöz
    Fh0 = 40;  % Yeni ceket sıvı akış hızı [m³/dak]
    y0_2 = y1(end,:); 
    [t2, y2] = ode15s(@CSTR_ODEs, t_final, y0_2);

    % çözümleri birleştir
    t = [t1; t2];
    y = [y1; y2];

    % Alt grafiklerle figure oluştur
    figure('Position', [100 100 800 600]);

    % Plot konsantrasyonu
    subplot(2,1,1);
    plot(t, y(:,1), 'LineWidth', 2);
    grid on;
    xlabel('Time (min)', 'FontSize', 12);
    ylabel('Concentration C_A (kmol/m^3)', 'FontSize', 12);
    title('Concentration Profile', 'FontSize', 14);
    % Add vertical line at flow rate change
    hold on;
    plot([20 20], [2 5.5], '--k', 'LineWidth', 1.5);
    text(21, 4, 'Flow rate change', 'FontSize', 10);
    xlim([0 40]);
    ylim([2 5.5]);

    % Plot sıcaklığı
    subplot(2,1,2);
    plot(t, y(:,2), 'LineWidth', 2);
    grid on;
    xlabel('Time (min)', 'FontSize', 12);
    ylabel('Temperature T (K)', 'FontSize', 12);
    title('Temperature Profile', 'FontSize', 14);
    % Akış hızı değişiminde dikey çizgi ekle
    hold on;
    plot([20 20], [295 335], '--k', 'LineWidth', 1.5);
    text(21, 315, 'Flow rate change', 'FontSize', 10);
    xlim([0 40]);
    ylim([295 335]);
end

function dydt = CSTR_ODEs(t, y)
    
    global F0 V CA0 T0 Cp rho k0 E R dH Th0 Cph rhoh Fh0
    
    
    CA = y(1);
    T = y(2);
    
    % Tepkime hızını hesapla
    k = k0 * exp(-E/(R*T));
    r = k * CA;
    
    % Isı üretimini hesapla
    Q = 12 * Fh0^2 * (Th0 - T);
    
    % Kütle dengesi
    dCAdt = (F0/V)*(CA0 - CA) - r;
    
    % Enerji dengesi
    dTdt = (F0/V)*(T0 - T) + Q/(rho*V*Cp) + (-dH)*r/(rho*Cp);
    
    % Dönüş türevleri
    dydt = [dCAdt; dTdt];
end