datafunc/+TwoPhaseData/Michel.m

classdef Michel < TwoPhaseData.Interface

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

      Skubus = S.*S.*S;
      kw = model.tp_kw_factor * Skubus / 2;
    end

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

      Squadrat = S.*S;
      dkw = model.tp_kw_factor * 3/2 * Squadrat;
    end

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

      Skubus = (1-S).*(1-S).*(1-S);
      ko = model.tp_ko_factor * Skubus ./ 3;
    end

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

      Squadrat = (1-S) .* (1-S);

      dko = (-1) .* model.tp_ko_factor .* Squadrat;
      assert(~any(isnan(dko)));
    end

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

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

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

      dpc = model.tp_pb .* (-1/2)./sqrt(S.^3.-S.^4);
    end

    function ddpc = capillary_pressure_second_derivative(this, glob, S, model)
      ddpc = model.tp_pb .* (1/4).*(S.^3-S.^4).^(-3/2).*(3.*S.^2-4.*S.^3);
%      ddpc = model.tp_pb .* (1/2).*(4./S - 1./(S.^2-S.^3));
%      ddpc = model.tp_pb .* (-1/2).*(-1.5.*S.^(0.5) + 2.*S)./(S.^(1.5)-S.^2).^2;
    end

    function c = injection_concentration(this, model)

      c = model.tp_injection_c;
    end

    function s_under = lower_source(this, elemin, loc, grid)

      glob = local2global(grid, elemin, loc, this);
      X = glob(:,1);
      Y = glob(:,2);

%      if model.t < 0.5
        D3 = ((X-0.8).^2 + (Y-0.5).^2) <= 0.01;
        s_under = 30 * D3;
%      else
%        s_under = zeros(size(X));
%      end
    end

    function s_above = upper_source(this, elemin, loc, grid)

      glob = local2global(grid, elemin, loc, this);
      X  = glob(:,1);
      Y  = glob(:,2);

%      if model.t < 0.5
      D1 = ((X-0.5).^2 + (Y-0.8).^2) <= 0.01;
      D2 = ((X-0.2).^2 + (Y-0.2).^2) <= 0.01;
      s_above = 10 * D1 + 20 * D2;
%      else
%        s_above = zeros(size(X));
%      end
    end

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