rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
Michel.m
1 classdef Michel < TwoPhaseData.Interface
2 
3  methods
4  function kw = water_permeability(this, glob, S, model)
5 
6  Skubus = S.*S.*S;
7  kw = model.tp_kw_factor * Skubus / 2;
8  end
9 
10  function dkw = water_permeability_derivative(this, glob, S, model)
11 
12  Squadrat = S.*S;
13  dkw = model.tp_kw_factor * 3/2 * Squadrat;
14  end
15 
16  function ko = oil_permeability(this, glob, S, model)
17 
18  Skubus = (1-S).*(1-S).*(1-S);
19  ko = model.tp_ko_factor * Skubus ./ 3;
20  end
21 
22  function dko = oil_permeability_derivative(this, glob, S, model)
23 
24  Squadrat = (1-S) .* (1-S);
25 
26  dko = (-1) .* model.tp_ko_factor .* Squadrat;
27  assert(~any(isnan(dko)));
28  end
29 
30  function pc = capillary_pressure(this, glob, S, model)
31 
32  pc = model.tp_pb * sqrt((1-S)./S);
33  end
34 
35  function dpc = capillary_pressure_derivative(this, glob, S, model)
36 
37  dpc = model.tp_pb .* (-1/2)./sqrt(S.^3.-S.^4);
38  end
39 
40  function ddpc = capillary_pressure_second_derivative(this, glob, S, model)
41  ddpc = model.tp_pb .* (1/4).*(S.^3-S.^4).^(-3/2).*(3.*S.^2-4.*S.^3);
42 % ddpc = model.tp_pb .* (1/2).*(4./S - 1./(S.^2-S.^3));
43 % ddpc = model.tp_pb .* (-1/2).*(-1.5.*S.^(0.5) + 2.*S)./(S.^(1.5)-S.^2).^2;
44  end
45 
46  function c = injection_concentration(this, model)
47 
48  c = model.tp_injection_c;
49  end
50 
51  function s_under = lower_source(this, elemin, loc, grid)
52 
53  glob = local2global(grid, elemin, loc, this);
54  X = glob(:,1);
55  Y = glob(:,2);
56 
57 % if model.t < 0.5
58  D3 = ((X-0.8).^2 + (Y-0.5).^2) <= 0.01;
59  s_under = 30 * D3;
60 % else
61 % s_under = zeros(size(X));
62 % end
63  end
64 
65  function s_above = upper_source(this, elemin, loc, grid)
66 
67  glob = local2global(grid, elemin, loc, this);
68  X = glob(:,1);
69  Y = glob(:,2);
70 
71 % if model.t < 0.5
72  D1 = ((X-0.5).^2 + (Y-0.8).^2) <= 0.01;
73  D2 = ((X-0.2).^2 + (Y-0.2).^2) <= 0.01;
74  s_above = 10 * D1 + 20 * D2;
75 % else
76 % s_above = zeros(size(X));
77 % end
78  end
79 
80  function descr = default_descr(this, descr)
81  descr.tp_pb = -0.5;
82  descr.tp_kw_factor = 0.15;
83  descr.tp_ko_factor = 1;
84  descr.tp_injection_c = 0.80;
85  end
86  end
87 end