Matlab backup

clear;
clc;

% Parâmetros
Sro = 0.2;        % Saturação residual do óleo
Srw = 0.2;        % Saturação residual da água
phi = 1;          % Porosidade
V = 1;            % Velocidade
a = 0.5;          % Multiplicador de viscosidade
L = 3;            % Comprimento do domínio (0,L)
T = 1;            % Tempo final
n = 100;          % Número de pontos
max_iter = 1000;  % Máximo de iterações
tol = 1e-6;       % Tolerância

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

dt = 0.001;
t = 0:dt:T;  % Malha temporal

% Parâmetros das permeabilidades relativas (Brooks & Corey)
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 fracionária de fluxo da água
fw = @(Sw) Krw(Sw) ./ (Krw(Sw) + a * Kro(Sw));

% Derivadas das funções
dKrw = @(Sw) krwm * nw * ((Sw - Srw).^(nw - 1)) / (1 - Srw - Sro).^nw;
dKro = @(Sw) -krom * no * ((1 - Sw - Sro).^(no - 1)) / (1 - Srw - Sro).^no;

dfw = @(Sw) a * (dKrw(Sw) .* Kro(Sw) - Krw(Sw) .* dKro(Sw)) ./ (Krw(Sw) + a * Kro(Sw)).^2;

% Inicialização da matriz B com as condições iniciais
B = Srw * ones(n, 1);
%B(1) = 1 - Sro;

% Inicialização da matriz solução
S = B;

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

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

S_k1 = S;

% Método de Newton
for j = 2:length(t)
    S_k = S;
    for iter = 1:max_iter
        for i = 1:n
            if i == 1
                % Função fracionária de fluxo da água
                fw_i = fw(S_k(i));
                dw_i = dfw(S_k(i));

                % Função g(S) e sua derivada g'(S)
                g = @(Sw) fw(Sw) + Sw * (alpha/dt);
                dg = @(Sw) dfw(Sw) + alpha/dt;

                % Atualização da saturação usando a equação (11)
                S_prev = S_k(i);
                for m = 1:max_iter
                    S_new = max(0, min(1, S_prev + (0.5 - g(S_prev)) / dg(S_prev))); % 0.5 é o valor de c_i utilizado
                    if abs(S_new - S_prev) < tol
                        break;
                    end
                    S_prev = S_new;
                end

                S_k1(i) = S_new;
            else
                % Função fracionária de fluxo da água
                fw_i_1 = fw(S_k(i - 1));
                fwi = fw(S_k(i));
                dw_i_1 = dfw(S_k(i - 1));
                dwi = dfw(S_k(i));

                % Função g(S) e sua derivada g'(S)
                g = @(Sw) fw(Sw) + Sw * (alpha/dt);
                dg = @(Sw) dfw(Sw) + alpha/dt;

                % Atualização da saturação usando a equação (11)
                S_prev = S_k(i);
                for m = 1:max_iter
                    S_new = max(0, min(1, S_prev + (0.5 - g(S_prev)) / dg(S_prev))); % 0.5 é o valor de c_i utilizado
                    if abs(S_new - S_prev) < tol
                        break;
                    end
                    S_prev = S_new;
                end

                S_k1(i) = S_new;
            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

    S = S_k;

    if ismember(t(j), tempos)
        plot(x, S, 'DisplayName', ['t = ', num2str(t(j))]);
        hold on;
    end
end

xlabel('x');
ylabel('Sw');
title('Solução Numérica');
legend;
grid on;
hold off;