BrooksCorey.m

classdef BrooksCorey < TwoPhaseData.Interface

  methods
    function kw = water_permeability(this, glob, S, model)

      Squadrat = S.*S;
      kw = model.tp_kw_factor * Squadrat ./ (Squadrat + (1-S).^(0.5));
    end

    function dkw = water_permeability_derivative(this, glob, S, model)

      oben = S .* S;
      unten = oben + (1-S).^(0.5);

      oben_der = 2 * S;
      unten_der = oben_der - 0.5 * (1-S).^(-0.5);

      dkw = model.tp_kw_factor * (oben_der .* unten - oben .* unten_der) ./ (unten .* unten);
      assert(~any(isnan(dkw)));
    end

    function ko = oil_permeability(this, glob, S, model)

      Squadrat = (1-S).*(1-S);
      ko = model.tp_ko_factor * Squadrat ./ (Squadrat + S.^(0.5));
    end

    function dko = oil_permeability_derivative(this, glob, S, model)

      oben = (1-S) .* (1-S);
      unten = oben + S.^(0.5);

      oben_der = -2 + 2 * S;
      unten_der = oben_der + 0.5 * S.^(-0.5);

      dko = model.tp_ko_factor * (oben_der .* unten - oben .* unten_der) ./ (unten .* unten);
      assert(~any(isnan(dko)));
    end

    function pc = capillary_pressure(this, glob, S, model)

      pc = model.tp_pb * S.^(-1/model.tp_lambda);
    end

    function dpc = capillary_pressure_derivative(this, glob, S, model)

      dpc = model.tp_pb .* (-1/model.tp_lambda) .* S.^(-1/model.tp_lambda - 1);
    end

    function ddpc = capillary_pressure_second_derivative(this, glob, S, model)

      ddpc = model.tp_pb .* (-1/model.tp_lambda).*(-1/model.tp_lambda - 1) .* S.^(-1/model.tp_lambda - 2);
    end

    function c = injection_concentration(this, glob, model)

      c = 0;
    end

    function s_under = lower_source(this, glob, model)

      s_under = 0;
    end

    function s_above = upper_source(this, glob, model)

      s_above = 0;
    end

    function descr = default_descr(this, descr)
      descr.tp_pb = -0.5;
      descr.tp_kw_factor = 1;
      descr.tp_ko_factor = 1;
      descr.tp_lambda = 0.5;
    end
  end
end