-
Notifications
You must be signed in to change notification settings - Fork 17
/
InputFile.m
136 lines (107 loc) · 7.41 KB
/
InputFile.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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
% ---------------------------------------------------------------------
% Copyright (C) 2016 by the Thermaid authors
%
% This file is part of Thermaid.
%
% Thermaid is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Thermaid is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Thermaid. If not, see <http://www.gnu.org/licenses/>.
% ---------------------------------------------------------------------
%
% Authors: Gunnar Jansen, University of Neuchatel, 2016-2017
% Ivan Lunati, Rouven Kuenze, University of Lausanne, 2012
%% GRID PARAMETERS ---------------------------------------------------------------------------------%
udata.len = [9 9]; % Physical length of the domain in x and y direction [m]
udata.Nf = [301 301]; % Number of cells in x and y direction
udata.dx = udata.len./udata.Nf; % Cell length [m]
%% SIMULATION PARAMETER FOR TRANSPORT --------------------------------------------------------------%
udata.timeSim = 1; % Total simulation time [s]
udata.dt = 1; % Time step length [s]
udata.tol = 1.e-4; % Tolerance on pressure-heat loop [-]
udata.maxit = 100; % Maximum number of pressure-heat loops to converge
%% FRACTURE NETWORK
% Crossed fracture test case
frac_cross
if (udata.dxf < min(udata.dx))
error('dxf < dx')
end
%% INITIAL CONDITIONS--------------------------------------------------------------------------------%
udata.T0 = 10*ones(udata.Nf(1),udata.Nf(2)); % Initial matrix temperature [°C]
udata.T0f = 10*ones(udata.Nf_f,1); % Initial fracture temperature [°C]
udata.tmax = 10; % Maximum temperature [°C] for plotting
udata.p0 = zeros(udata.Nf(1),udata.Nf(2)); % Initial matrix pressure [Pa]
udata.p0f = zeros(udata.Nf_f,1); % Initial fracture pressure [Pa]
%% BC FLUID ----------------------------------------------------------------------------------------%
udata.ibcs = zeros(2*sum(udata.Nf),1); % Type 0:Neumann(N); 1:Dirichlet(D)
udata.Fix = zeros(2*sum(udata.Nf),1); % Value N [m2/s] (inflow>0); D [Pa]
udata.ibcs(1:udata.Nf(2)) = 1;
udata.ibcs(udata.Nf(2)+1:2*udata.Nf(2))=1;
udata.Fix(1:udata.Nf(2)) = 1;
udata.Fix(udata.Nf(2)+1:2*udata.Nf(2))=0;
%% BC TRANSPORT ------------------------------------------------------------------------------------%
udata.flagHeatTransport = 0; % If 1 - heat transport is used
udata.FixT = zeros(2*sum(udata.Nf),1); % Temperature at the boundary [°C]
%% SOURCE TERMS ------------------------------------------------------------------------------------%
Q = zeros(udata.Nf); % Source term [m2/s]; inflow positive
QT = zeros(udata.Nf); % Inflow temperature [°C]
%% GRAVITY---------------------------------------------------------------------------------%
udata.gravity = 0; % Gravity acceleration in y-direction [m/s2]
%% PERMEABILITY ------------------------------------------------------------------------------------%
global k_ratio
a = 1/250;
udata.K = ones(udata.Nf(1),udata.Nf(2))*1; % Permeability field [m2]
udata.K_f = ones(udata.Nf_f,1)*k_ratio; % Fracture permeability field [m2]
%% Fracture aperture
udata.b0 = a*ones(udata.Nf_f,1); % Fracture aperture field [m]
%% Porosity ----------------------------------------------------------------------------------------%
udata.phi = ones(udata.Nf(1),udata.Nf(2))*0.3; % Porosity field [-]
udata.phi_f = ones(udata.Nf_f,1)*0.3; % Fracture porosity field [-]
%% Fluid density ------------------------------------------------------------------------%
udata.const_density = 1; % If 1 - constant fluid density is used
if(udata.const_density)
udata.density_l = 1000*ones(udata.Nf); % Density of the fluid [kg/m3]
udata.density_lf = 1000*ones(udata.Nf_f,1); % Density of the fluid [kg/m3]
end
%% Fluid viscosity ------------------------------------------------------------------------%
udata.const_viscosity = 1; % If 1 - constant fluid viscosity is used
if(udata.const_viscosity)
udata.viscosity = 1e-3*ones(udata.Nf); % Viscosity of the fluid [Pa*s]
udata.viscosity_f = 1e-3*ones(udata.Nf_f,1); % Viscosity of the fluid [Pa*s]
end
%% ROCK DENSITY ------------------------------------------------------------------------%
udata.density_s = 2000*ones(udata.Nf); % Density of the rock [kg/m3]
udata.density_sf = 2000*ones(udata.Nf_f,1); % Density of the rock [kg/m3]
%% MECHANIC PROPERTIES
udata.flagIncompressible = 1; % If 1 - incompressible fluids are used
udata.compressibility_l = 5e-1000; % Compressibility of the fluid [1/Pa]
udata.compressibility_s = 5e-1000; % Compressibility of the rock [1/Pa]
udata.shear_modulus = 29e9; % Shear modulus of the rock [Pa]
udata.poisson_ratio = 0.25; % Poisson ratio of the rock [-]
udata.friction_coeff = 0.6*ones(udata.Nf_f,1); % Friction coefficient (shear failure) [-]
udata.K_enh = 1e2; % Permeability enhancement factor [-]
udata.therm_exp_coeff = 7.9e-6; % Thermal expansion coefficient of the rock matrix [-]
udata.flagFracStability = 0; % If 0 the remaining parameters in this section are not used
udata.sigma_1 = 20e6*ones(udata.Nf_f,1); % Maximum principal stress (S1) [Pa]
udata.sigma_2 = 14e6*ones(udata.Nf_f,1); % Intermeadiate principal stress (S2) [Pa]
udata.sigma_3 = 12e6*ones(udata.Nf_f,1); % Minimum principal stress (S3) [Pa]
% TR - Trend (from N) / PL - Plunge (from horizontal)
udata.stress_trend = [0 90]; % [TR(S1) TR(S3)] [°]
udata.stress_plunge = 90; % PL(S1) [°]
udata.frac_az = 90*ones(size(udata.frac_angle)); % Fracture dip [°]
udata.frac_dip = min(abs(udata.frac_angle),180-abs(udata.frac_angle)); % Fracture azimuth (from N) [°]
udata.use_thermal_stress = 0; % Boolean for thermal stress in calculations. If 0 the contribution of thermal stress is not calculated.
%% THERMAL DIFFUSION ------------------------------------------------------------------------%
udata.lambda_l = 0; % Thermal conductivity of the fluid [W/(m*K)]
udata.lambda_s = 0; % Thermal conductivity of the rock [W/(m*K)]
udata.cp_l = 0; % Specific heat capacity of the fluid [J/(kg*K)]
udata.cp_s = 0; % Specific heat capacity of the rock [J/(kg*K)]
udata.ibcD = zeros(2*sum(udata.Nf),1); % 1 -> Diffusion on boundary cells