How to Set Zero Transmissibility in MRST for Flow Isolation?

[Calvin Scott]

Dear all,

I'm currently working on fitting a model to experimental data. To accurately replicate the experimental setup, I need to set the transmissibility across a specific face within the reservoir (perpendicular to the horizontal well direction) to zero (or a very small value). The goal is to effectively divide the reservoir into two isolated regions with no fluid exchange between them.

However, I'm unsure how to directly set or modify the transmissibility in my model and incorporate this change into the simulation workflow. I would greatly appreciate any guidance or advice on how to achieve this.

Thank you.

Below is my model script:

%%%%%%%%%%%%%%%%%%%% mrstModule add ad-core ad-blackoil ad-props test-suite deckformat wellpaths nwm hfm mrst-gui streamlines

% 使用井拓扑建立多段井模型

%% 基质模型

fn = fullfile('tupdate_Trans.DATA');

deck = readEclipseDeck(fn);

deck = convertDeckUnits(deck); % 包含岩石属性，流体属性等源Eclipse data中的各个元素

fluid = initDeckADIFluid(deck);

% fluid = assignRelPerm(fluid); % assignRelPerm用于油气水三相系统

%% SETUP GRID

x = 0:20:1440;

y = unique([0:20:1000]); % 490:0.1:510]);

% y = [0:20:480 480:1:520 520:20:1000];

% y = unique(y);

z = 4900:20:4960;

G = computeGeometry(tensorGrid(x, y, z)); % 原始网格储层结构——未添加实际深度

nx = G.cartDims(1);

ny = G.cartDims(2);

nz = G.cartDims(3);

%% set up permeability based on I-indices

%% 赋值部分网格

% 逻辑索引模型网格G在X,Y,Z轴的部分值

[I,J,K] = gridLogicalIndices(G);

I0 = I <= 4 ;

I1 = I >= 5 & I <= 12 ;

I2 = I >= 13 & I <= 20;

I3 = I >= 21 & I <= 28;

I4 = I >= 29 & I <= 36;

I5 = I >= 37 & I <= 44;

I6 = I >= 45 & I <= 52;

I7 = I >= 53 & I <= 60;

I8 = I >= 61 & I <= 68;

I00 = I >= 69 & I <= 72;

px = ones(G.cells.num,1);

px(I0) = 1.4 *milli*darcy;

px(I1) = 1.4 *milli*darcy;

px(I2) = 1.4 *milli*darcy;

px(I3) = 1.4 *milli*darcy;

px(I4) = 1.4 *milli*darcy;

px(I5) = 1.4 *milli*darcy;

px(I6) = 7.28 *milli*darcy;

px(I7) = 0.17 *milli*darcy;

px(I8) = 0.17 *milli*darcy;

px(I00)= 0.17 *milli*darcy;

% 非均质孔隙度分布情况

poro = ones(G.cells.num,1);

poro(I0) = 0.09 ;

poro(I1) = 0.09 ;

poro(I2) = 0.09 ;

poro(I3) = 0.09 ;

poro(I4) = 0.09 ;

poro(I5) = 0.09 ;

poro(I6) = 0.1 ;

poro(I7) = 0.08 ;

poro(I8) = 0.08 ;

poro(I00) = 0.08 ;

c = (nx*ny*(1+1/2)+5):(nx*ny*(1+1/2)+68);

% 定义目标渗透率及对应的 skin

perm_value = [1.4, 7.28, 0.17];

skins_val = [4, -2.5, 2.5];

skin = autoSkinFlexible(c, px, 'value', perm_value, skins_val);

% 清理掉参数c等临时参数

clear c perm_value skins_val;

%% 再方案设置位置修改文件名，方便模块化修改方案

excelFilename = 'BAC-512.xlsx';

%% --------------------------------------

G = computeGeometry(G);

perm = px;

% 将rock属性置于地质结构 G 中

G.rock = makeRock(G,px,poro);

%% Define initial composition and pressure 定义初始组分和初始压力

% 定义初始油藏压力

gravity on; % tell MRST to turn on gravity

g = gravity; % get the gravity vector-注意重力加速度是三维向量

rhow = 1000; % density of brine corresponding to 94 degrees C and 300 bar. 94℃和300bpa压力下的盐水密度

initState.pressure = rhow*g(3)*G.cells.centroids(:,3); % 计算静水压力分布

% My reservoir only contains gas and water.

% 定义流体组分-包含岩石属性 初始化储层 水 油 气

s0 = [0.3, 0, 0.7];

state0 = initResSol(G, initState.pressure, s0);

%% Find fracture cell faces at domain boundary

% boundfaces=findfracboundaryfaces(G,tol);

%% 设置压力边界值 fp (fpz表示纵向边界压力值) —— 其下都是采用网格索引

% 从一组网格单元中提取边界面

% f 表示提取的边界面，包括这些面在G数组中的索引位置————定位

f = boundaryFaces(G);

fn = G.faces.normals(f,:);

fnz = G.faces.normals(f,3);

f = f(abs(G.faces.normals(f,3))<eps);

fpz = G.faces.centroids(f,3)*rhow*norm(gravity);

bc = [];

bc = addBC(bc,f,'pressure',fpz,'sat',[0 0 1]);

%% relative permeability and capillary pressure

% 三相中气水系统的相对渗透率表 Sg Krg Krw

SWGF = {

[[ 0.08 0 0.0665 ] ;

[ 0.302 0.005 0.012 ] ;

[ 0.34 0.01 0.006 ] ;

[ 0.4 0.017 0.003 ] ;

[ 0.415 0.02 0.0023 ] ;

[ 0.433 0.024 0.0015 ] ;

[ 0.47 0.033 0.0012 ] ;

[ 0.55 0.07 0 ] ;

[ 0.73 0.15 0 ] ;

]};

% 三相中油水系统的相对渗透率表 Sw Krw Kro

SWOF = {

[[ 0.08 0 0.0665 ] ;

[ 0.302 0.005 0.012 ] ;

[ 0.34 0.01 0.006 ] ;

[ 0.4 0.017 0.003 ] ;

[ 0.415 0.02 0.0023 ] ;

[ 0.433 0.024 0.0015 ] ;

[ 0.47 0.033 0.0012 ] ;

[ 0.55 0.07 0 ] ;

[ 0.73 0.15 0 ] ;

]};

%

krg = @(sg) interpTable(SWGF{1}(:,1),SWGF{1}(:,2),sg) ;

krw = @(sw) interpTable(SWGF{1}(:,1),SWGF{1}(:,3),1-sw) ;

kro = @(so) interpTable(SWOF{1}(:,1),SWOF{1}(:,3),1-so) ;

krow = @(so) interpTable(SWOF{1}(:,1),SWOF{1}(:,3),1-so) ;

krog = @(sg) interpTable(SWGF{1}(:,1),SWGF{1}(:,3),1-sg) ;

% 修改相渗端点值

fluid.krW = krw;

fluid.krO = kro;

fluid.krG = krg;

fluid.krOW = krow;

fluid.krOG = krog;

% 绘制相渗曲线

figure;

plot(SWGF{1}(:,1),...

SWGF{1}(:,[2,3]),'*-',...

'LineWidth', 2, 'MarkerSize', 5)

legend({'krg','krw'}, 'Location', 'Best')

xlabel('sg'), title('Relative Permeability')

%-----------------------------------------------------

%============================================================================================

% Define three-phase compressible flow model

model = ThreePhaseBlackOilModel(G,G.rock, fluid, 'disgas', false, 'vapoil', false);

c = (nx*ny*(1+1/2)+5):(nx*ny*(1+1/2)+68);

nc= numel(c);

prodS = addWell([],G,G.rock,c,'type','bhp','Val', 433*barsa,'Radius',.06*meter,'Dir','x','Skin',skin,'name','prod', 'compi', [0, 0, 1]);

prodS.lims.bhp = 100*barsa;

% Define additional properties needed for ms-well

% We have 8 node-to-node segments

topo = [1 2 3 4 5 6 7 8

2 3 4 5 6 7 8 9]' ;

cell2node = sparse(rldecode((2:9)',8*ones(8,1)),(1:64)',1,9,64); % 可在permTensor.m函数43行看到rldecode用法

spy(cell2node);

% Segment lengths/diameters and node depths/volumes

segments = 8;

lengths = 20*8*ones(8,1);

diam = .1*ones(8,1);

depths = G.cells.centroids(c([1 1:segments]),3) ;

vols = ones(9,1);

% Convert to ms-well

prodMS = convert2MSWell(prodS, 'cell2node', cell2node, 'topo', topo, 'G', G, 'vol', vols, ...

'nodeDepth', depths, 'segLength', lengths, 'segDiam', diam);

% Finally, we set flow model for each segment:

% segments 1-9: wellbore friction model

% no pressure relief valve

wbix = deal(1:8) ;

roughness = 1e-3; % discharge 表示流量系数，可以用此来表示阀门

% Set up flow model as a function of velocity, density and viscosity

prodMS.segments.flowModel = @(v, rho, mu)...

[wellBoreFriction(v(wbix), rho(wbix), mu(wbix), prodMS.segments.diam(wbix), ...

prodMS.segments.length(wbix), roughness, 'massRate')];

% 使用物理归一化残差进行非线性收敛计算。

model.useCNVConvergence = true;

% Setting up the non-linear solver.

nonlinearsolver = NonLinearSolver();

nonlinearsolver.useRelaxation = true;

% Simple uniform schedule with initial rampup

totTime = 900*day;

dt = rampupTimesteps(totTime, 30*day);

% 额外添加了边界条件的schedule

schedule = simpleSchedule(dt, 'W', prodMS,'bc',bc);

% 未额外添加边界条件的schedule

% schedule = simpleSchedule(dt, 'W', prodMS);

[ws, states,reports] = simulateScheduleAD(state0, model, schedule, 'NonLinearSolver', nonlinearsolver);

[Olav Møyner]

Hi,

You can modify model.operators.T. You can look up what pairs of cells are connected to connection face by model.operators.N.

Bset regards,

Olav

---

From: sinte...@googlegroups.com <sinte...@googlegroups.com> on behalf of Calvin Scott <kungcli...@gmail.com>
Sent: Friday, May 9, 2025 19:58
To: MRST-users: The Matlab Reservoir Simulation Toolbox User Group <sinte...@googlegroups.com>
Subject: [MRST Users] How to Set Zero Transmissibility in MRST for Flow Isolation?

 

--
You received this message because you are subscribed to the Google Groups "MRST-users: The Matlab Reservoir Simulation Toolbox User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sintef-mrst...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/sintef-mrst/1bee2585-b6d4-424c-a004-191e4d5953c8n%40googlegroups.com.