How to build field-scale fractured geothermal model with MRST?

150 views

Skip to first unread message

Chen Guodong

unread,

Sep 2, 2022, 3:43:29 AM9/2/22

to sinte...@googlegroups.com, mr...@sintef.no, olav....@sintef.no, bard.sk...@sintef.no

Dear MRST team,

Thanks for develop an outstanding open-source simulation platform. Recently I was trying to construct 3D field-scale fractured geothermal model with MRST. The fractured flow simulation was achieved successfully. However, after combining hfm and geothermal module, the built model cannot run successfully. Is MRST able to achieve 3D field-scale fractured geothermal simulation and do you have such example?

%% LOAD MRST MODULES

clc,clear

mrstModule add hfm ad-blackoil ad-core ad-props compositional geothermal mrst-gui

% path1 = 'D:\PhD_research\Multifractal_characterization\mrst2021b';

% addpath(genpath(path1));

path2 = 'D:\PhD_research\Multifractal_characterization\shale_main';

addpath(genpath(path2));

%% FRACTURE NETWORK STATISTICAL PARAMETERS

physdim=[1000 500 100]; % domain size

% Fracture sets 1 (red), 2(green), 3(yellow)

fracinput1=struct('vertices',10,...

'P32',0.0005*(1/meter),...

'normal',struct('direction',[1 0 0],'K',10),...

'size',struct('minsize',10*meter,'maxsize',50*meter,'exponent',1.5),...

'perm',1e7*milli*darcy,...

'poro',1,...

'aperture',0.001*meter,...

'circle',true);

fracinput2=fracinput1; fracinput2.normal.direction=[0 0 1];

fracinput3=fracinput1; fracinput3.normal.direction=[0 1 0];

% Exclusion zone is cylindrical with radius being (1+exclzonemult) times

% the fracture radius; height is exclzonemult times the fracture radius.

exclzonemult=0.001;

%% DFN GENERATION

tol=10^-5; % tolerance for comparison of doubles

fracplanes=[]; % empty list of fracture objects

% set 1

[fracplanes,~]=DFNgenerator(fracplanes,fracinput1,physdim,...

exclzonemult,tol);

% set 2

[fracplanes,~]=DFNgenerator(fracplanes,fracinput2,physdim,...

exclzonemult,tol);

% set 3

[fracplanes,~]=DFNgenerator(fracplanes,fracinput3,physdim,...

exclzonemult,tol);

%% SETUP GRID

celldim=[50 25 10];

G = cartGrid(celldim, physdim);

G = computeGeometry(G);

%% Make rock structure

km=10*milli*darcy;

G.rock=makeRock(G,km,0.25);

% figure; plotGrid(G,'facealpha',0);

% view(15,20);

% Add thermal props

G.rock = addThermalRockProps(G.rock, 'lambdaR', 2, 'rhoR', 2700, 'CpR', 1000);

%% PLOT DFN

colourchoice=['r','g','y'];

figure;

hold on;

plotGrid(cartGrid([1 1 1],physdim),'facealpha',0);

axis equal tight

view(45,30);

for i=1:3

index = find(vertcat(fracplanes.SetID)==i);

C=colourchoice(i);

for j=index'

X=fracplanes(j).points(:,1);

Y=fracplanes(j).points(:,2);

Z=fracplanes(j).points(:,3);

fill3(X,Y,Z,C);

end

xlabel('x','Interpreter','latex')

ylabel('y','Interpreter','latex')

zlabel('z','Interpreter','latex')

end

%% EDFM PRE-PROCESSING

tol=1e-6;

[G,fracplanes]=EDFMgrid(G,fracplanes,'Tolerance',tol);

G=fracturematrixNNC3D(G,tol);

[G,fracplanes]=fracturefractureNNCs3D(G,fracplanes,tol);

%% SETUP FLUID MODEL WITH WATER PROPERTIES

pRef = 100*barsa;

fluid = initSimpleADIFluid('phases','W',...

'mu' , 1* centi*poise , ...

'rho', 1000* kilogram/meter^3, ...

'pRef', pRef, ...

'c', 1e-5/barsa);

fluid = addThermalFluidProps(fluid, 'Cp' , 4.2e3, ...

'lambdaF', 0.6 , ...

'useEOS' , true );

%% Inititate state transient flow

state = initResSol(G, pRef); % Pressure and saturation

state.T = ones(G.cells.num,1).*(273.15+10); % Temperature

%% Define the problem - GeothermalWaterModel for pure water problem

% Define the gravity vector for a 3D x,y,z case using MRST convention

gravity reset on

% gravity([0 0 -9.81]); % gravity is downward

model = GeothermalModel(G, G.rock, fluid, [], 'extraStateOutput', true, 'outputFluxes', true);

% model.FlowDiscretization = setTimeDiscretization(model, 'fully-implicit'); % FIM

% operators = model.operators;

% model.operators = setupEDFMOperatorsTPFA(G, G.rock, tol);

% model.operators.AccDiv = operators.AccDiv;

% @(acc,flux)acc+C'*flux

% model.stepFunctionIsLinear=true;

% Define limits of temperature validity for EOS

K0 = 273.15*Kelvin;

TMin = K0; % Minimum temperature

TMax = K0 + 275*Kelvin; % Maximum temperature

pMax = 200e6*Pascal; % Maximum pressure

model.minimumTemperature = TMin;

model.maximumTemperature = TMax;

model.maximumPressure = pMax;

%% ADD INJECTOR

% totTime = 5*year;

totTime = 3000*day;

tpv = sum(model.operators.pv);

wellRadius = 0.1;

[nx, ny, nz] = deal(G.cartDims(1), G.cartDims(2), G.cartDims(3));

cellinj1 = nx/2+1:nx*ny: nz*nx*ny;

W = addWell([], G, G.rock, cellinj1, 'Type', 'rate', ...

'Val', tpv/totTime, 'Radius', wellRadius, ...

'Comp_i', [1, 0, 0], 'Name', 'Injector1');

cellinj2 = nx*ny/2-nx/2+1:nx*ny: nz*nx*ny;

W = addWell(W, G, G.rock, cellinj2, 'Type', 'rate', ...

'Val', tpv/totTime, 'Radius', wellRadius, ...

'Comp_i', [1, 0, 0], 'Name', 'Injector2');

cellinj3 = nx*ny/2+nx/2: nx*ny : nz*nx*ny;

W = addWell(W, G, G.rock, cellinj3, 'Type', 'rate', ...

'Val', tpv/totTime, 'Radius', wellRadius, ...

'Comp_i', [1, 0, 0], 'Name', 'Injector3');

cellinj4 = nx*ny-nx/2+1 : nx*ny : nz*nx*ny;

W = addWell(W, G, G.rock, cellinj4, 'Type', 'rate', ...

'Val', tpv/totTime, 'Radius', wellRadius, ...

'Comp_i', [1, 0, 0], 'Name', 'Injector4');

%% ADD PRODUCER

cellprod1 = 1 : nx*ny : nz*nx*ny;

W = addWell(W, G, G.rock, cellprod1, 'Type', 'bhp', ...