Matlab Backup 2

clear;
clc;

% Parâmetros globais
Sro = 0;      % Saturação residual do óleo
Srw = 0;     % Saturação residual da água
phi = 1;          % Porosidade
V = 1;            % Velocidade
a = 0.5;         % Multiplicador de viscosidade
L = 3;            % Comprimento do domínio
T = 1;            % Tempo percorrido

% Função Krw e Kro
nw = 2;
no = 2;
krwm = 1;
krom = 1;
Krw = @(Sw) krwm .* ((Sw - Srw) ./ (1 - Srw - Sro)).^nw;
Kro = @(Sw) krom .* ((1 - Sw - Sro) ./ (1 - Srw - Sro)).^no;

% Função fw
fw = @(Sw) Krw(Sw) ./ (Krw(Sw) + a * Kro(Sw));

% Derivadas das funções
dfw = @(Sw) (krwm*((Srw - Sw)/(Sro + Srw - 1))^nw*((krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 1))/(Sro + Srw - 1) - (a*krom*no*((Sro + Sw - 1)/(Sro + Srw - 1))^(no - 1))/(Sro + Srw - 1)))/(krwm*((Srw - Sw)/(Sro + Srw - 1))^nw + a*krom*((Sro + Sw - 1)/(Sro + Srw - 1))^no)^2 - (krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 1))/((krwm*((Srw - Sw)/(Sro + Srw - 1))^nw + a*krom*((Sro + Sw - 1)/(Sro + Srw - 1))^no)*(Sro + Srw - 1));
d2fw = @(Sw) (2*krwm*((Srw - Sw)/(Sro + Srw - 1))^nw*((krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 1))/(Sro + Srw - 1) - (a*krom*no*((Sro + Sw - 1)/(Sro + Srw - 1))^(no - 1))/(Sro + Srw - 1))^2)/(krwm*((Srw - Sw)/(Sro + Srw - 1))^nw + a*krom*((Sro + Sw - 1)/(Sro + Srw - 1))^no)^3 - (krwm*((Srw - Sw)/(Sro + Srw - 1))^nw*((krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 2)*(nw - 1))/(Sro + Srw - 1)^2 + (a*krom*no*((Sro + Sw - 1)/(Sro + Srw - 1))^(no - 2)*(no - 1))/(Sro + Srw - 1)^2))/(krwm*((Srw - Sw)/(Sro + Srw - 1))^nw + a*krom*((Sro + Sw - 1)/(Sro + Srw - 1))^no)^2 + (krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 2)*(nw - 1))/((krwm*((Srw - Sw)/(Sro + Srw - 1))^nw + a*krom*((Sro + Sw - 1)/(Sro + Srw - 1))^no)*(Sro + Srw - 1)^2) - (2*krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 1)*((krwm*nw*((Srw - Sw)/(Sro + Srw - 1))^(nw - 1))/(Sro + Srw - 1) - (a*krom*no*((Sro + Sw - 1)/(Sro + Srw - 1))^(no - 1))/(Sro + Srw - 1)))/((krwm*((Srw - Sw)/(Sro + Srw - 1))^nw + a*krom*((Sro + Sw - 1)/(Sro + Srw - 1))^no)^2*(Sro + Srw - 1));

% % Ponto de inflexão
Sc = fsolve(d2fw, 0.5);

% Parâmetros para a solução numérica
n = 200;         % Número de pontos
max_iter = 100;  % Máximo de iterações
tol = 1e-6;      % Tolerância

dx = L / (n - 1);
x_num = linspace(dx/2, L-dx/2, n-1);  % Malha espacial centralizada

dt = 100 / 1000;
t = 0:dt:T;  % Malha temporal

% Inicialização das soluções anteriores com as condições iniciais
S_ant = Srw * ones(n-1, 1);

% Parâmetro da solução numérica
alpha = dx / V;

% Tempos que serão plotados
tempos = [0.5, 1];

% Método do Ponto Fixo Modificado
S_k1 = zeros(n-1,1);
for j = 2:length(t)
    S_k = Sc*ones(n-1,1);

    for iter = 1:max_iter
        for i = 1:n-1
            if i == 1
                beta = fw(1-Sro) / alpha;

                g = fw(S_k(i)) + (alpha / dt) * S_k(i);  % Função g
                dg = dfw(S_k(i)) + alpha / dt;  % Função g'

                gama = alpha * (beta + (S_ant(i)) / dt);

                % Atualização da saturação
                S_new = S_k(i) + (gama - g) / dg;
                S_k1(i) = max(0, min(1, S_new));
            else
                beta = (fw(S_k(i-1)) + dfw(S_k(i-1)) * (S_k1(i-1) - S_k(i-1))) / alpha;

                g = fw(S_k(i)) + (alpha / dt) * S_k(i);
                dg = dfw(S_k(i)) + alpha / dt;

                gama = alpha * (beta + (S_ant(i)) / dt);

                % Atualização da saturação
                S_new = S_k(i) + (gama - g) / dg;
                S_k1(i) = max(0, min(1, S_new));
            end

            % Passo que garante a incondicionalidade da convergência
            if d2fw(S_k1(i)) * d2fw(S_k(i)) < 0
                %S_k1(i) = (S_k1(i) + S_k(i))/2;
                S_k1(i) = Sc;
            end
        end

        % Verificação do critério de convergência
        if norm((S_k1 - S_k) / S_k1) < tol
            break;
        end

        S_k = S_k1;
    end

    % Atualização da solução anterior
    S_ant = S_k1;

    % Plot das soluções em tempos especificados
    if ismember(t(j), tempos)
        plot(x_num, S_ant, 'DisplayName', sprintf('Numérica t = %.2f', t(j)));
        hold on;
    end
end

% Malha espacial da analítica
dx_ana = 0.001;
x_ana = 0:dx_ana:L;

% Malha temporal da analítica
dt_ana = 0.1 / 1000;
t_ana = 0:dt_ana:T;

% Função Krw e Kro para a solução analítica
Krw_ana = @(Sw) krwm * ((Sw - Srw) / (1 - Srw - Sro)).^nw;
Kro_ana = @(Sw) krom * ((1 - Sw - Sro) / (1 - Srw - Sro)).^no;
fw_ana = @(Sw) Krw_ana(Sw) ./ (Krw_ana(Sw) + a * Kro_ana(Sw));
dKrw_ana = @(Sw) krwm * nw * ((Sw - Srw).^(nw - 1)) / (1 - Srw - Sro).^nw;
dKro_ana = @(Sw) -krom * no * ((1 - Sw - Sro).^(no - 1)) / (1 - Srw - Sro).^no;
dfw_ana = @(Sw) a * (dKrw_ana(Sw) .* Kro_ana(Sw) - Krw_ana(Sw) .* dKro_ana(Sw)) ./ (Krw_ana(Sw) + a * Kro_ana(Sw)).^2;

% Função para encontrar a raiz
func = @(Sw) dfw_ana(Sw) .* (Sw - Srw) - fw_ana(Sw) + fw_ana(Srw);

% Encontra a raiz de func(Sw)
Chute_inicial = 0.7;
Sw_choque = fsolve(func, Chute_inicial);

% Loop principal para a solução analítica
for t_val = tempos
    Sw_rf = linspace(1 - Sro, Sw_choque, 1000);
    x_rf = zeros(1, length(Sw_rf));

    x_choque = (V * t_val / phi) * dfw_ana(Sw_choque);
    for k = 1:length(Sw_rf)
        x_rf(k) = (V * t_val / phi) * dfw_ana(Sw_rf(k));
    end

    sw = [Sw_rf, Srw, Srw];
    xw = [x_rf, x_choque, L];
    plot(xw, sw, 'DisplayName', sprintf('Analítica t = %.2f', t_val));
    hold on;
end

xlabel('x');
ylabel('Sw');
title('Comparação da Solução Numérica e Analítica de Buckley-Leverett');
legend('show');
grid on;
hold off;