-
Notifications
You must be signed in to change notification settings - Fork 0
/
getSmoothedFracFlow.m
61 lines (58 loc) · 2.12 KB
/
getSmoothedFracFlow.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
%%%Copyright 2018 SINTEF AS
function varargout = getSmoothedFracFlow(model, p0, ph, ix)
if nargin < 3
ph = 'W';
end
if nargin < 4
ix = 1;
end
ap = model.getActivePhases;
switch ph
case 'W'
f = model.fluid;
p_av = mean(p0);
if isfield(model.fluid, 'sWcon')
sMin = model.fluid.sWcon;
else
sMin = 0;
end
if numel(sMin) ~= 1
sMin = sMin(ix);
end
sw = (sMin:.01:1)';
if ~ap(3) % no gas
[kw, ko] = model.relPermWO(sw, 1-sw, f);
[lw, lo] = deal(kw/f.muW(p_av), ko/f.muO(p_av));
else % assume we are interested in oil/water flow ...
[kw, ko] = model.relPermWO(sw, 1-sw, f, 'cellInx', repmat(ix, [numel(sw), 1]));
rs_sat = f.rsSat(p_av, 'cellInx', ix);
[lw, lo] = deal(kw/f.muW(p_av, 'cellInx', ix), ko/f.muO(p_av, rs_sat, true, 'cellInx', ix));
end
f = lw./(lw+lo);
% fit smooth function of the form a*sc(s)^n/(a*sc(s)^n - (1-sc(s))^n)
sc = @(s)max(0, s - sMin)./(1-sMin);
% coefficent function for given n (minimize misfit by least square)
coeff = @(n)((f-1).*sc(sw).^n)\(-f.*(1-sc(sw)).^n);
% smoothed fracflow of samples
ff = @(n)(coeff(n)*sc(sw).^n)./(coeff(n)*sc(sw).^n + (1-sc(sw)).^n);
% find best exponent 1<= n <= 5
fit = @(n)norm(ff(n)-f);
[n, ~, flag] = fminbnd(fit, 1, 5);
if flag > 0
a = coeff(n);
varargout{1} = @(s)(a*sc(s).^n)./(a*sc(s).^n + (max(0,1-sc(s))).^n);
if nargout > 1
num = @(s)a*sc(s).^n;
dnum = @(s)(a*n/(1-sMin))*sc(s).^(n-1);
den = @(s)(a*sc(s).^n + (max(1-sc(s), 0)).^n);
dden = @(s)(a*n/(1-sMin))*sc(s).^(n-1) - (n/(1-sMin))*(max(0, 1-sc(s))).^(n-1);
varargout{2} = @(s)(dnum(s).*den(s)-num(s).*dden(s))./(den(s).^2);
end
else
error('Unable to produce frac-flow curves')
%varargout = getFracFlow(model, state, ph);
end
otherwise
error('Not implemented')
end
end