Stress Maps application added

This commit is contained in:
Joanna Kocot 2017-03-31 13:00:31 +02:00
commit 5ecc6387c0
2 changed files with 597 additions and 0 deletions

4
.gitignore vendored Normal file
View File

@ -0,0 +1,4 @@
/*
/*/
!.gitignore
!/src/

View File

@ -0,0 +1,593 @@
function stress_maps_service(root_folder, resolution, weightings, keyword, nparam, mc, stresstype, ldvalue, plot_fracs)
% STRESS_MAPS_SERVICE Plot stress patterns associated with
% the hydraulic fracturing of a shale gas well. The data is based on
% the physical and geological conditions at the Preese Hall Episode site.
% STRESS_MAPS_SERVICE(ROOT_FOLDER, RESOLUTION, WEIGHTINGS, KEYWORD, NPARAM, MC, STRESSTYPE, LDVALUE)
% plots a stress map at a high ([20]m) or low ([50]m) RESOLUTION for the STRESSTYPE
% ([N]ormal, [S]hear or [C]oulomb) using a given source WEIGHTING ratio. The parameter
% to be plotted is defined in the variable KEYWORD (ISIP, FlowRate, PumpTime) and the value
% of this in NPARAM. MC is the monte carlo iteration (1-50) and ROOT_FOLDER the location
% of all the mat files. The lateral distance is calculated using the
% LDVALUE threshold
% Input options for nparam:
% dP = [0.426 0.676 0.926 1.176 1.676 1.926 2.176 2.426 2.926 3.426 3.676 4.426 4.926 5.176 6.676 8.176 9.676 11.176 12.676 14.176 15.676 17.176 18.676 20.176 21.676 23.176 24.676 26.176 27.676 29.176]
% FlowRate = [0.0983 0.1099 0.1215 0.1331 0.1446 0.1562 0.1678 0.1793 0.1909]
% PumpTime = [30 60 90 120 150 180 210 240 270 300 330 360]
% version 1.1 20/12/2016
% Copyright 2016 Rachel Westwood, Keele University
% This work is licensed under a Creative Commons Attribution-ShareAlike
% 4.0 International License. For full details of the license see
% https://creativecommons.org/licenses/by-sa/4.0/legalcode
% Version History
% 1.1 - corrected colorbar issue with 20m resolution and filename issue
% with DP
% check the weightings to ensure they are correct
if length(weightings) ~= 3
error('source weighting combination must contain 3 elements');
end
if sum(weightings) ~= 1
error('source weighting combination must add up to 1');
end
if ldvalue < 0.001 || ldvalue > 1
error('Lateral respect distance stress threshold should be between 0.001 and 1 MPa');
end
% check that nparam is in the list
dP = round(([37.25 37.5 37.75 38 38.5 38.75 39 39.25 39.75 40.25 40.5 41.25 41.75 42:1.5:66] - 36.824)*1000)/1000; % dP in MPa
flow = round((8500:1000:16500) * 1.157e-5 * 10000)/10000; % flow in m3/s
pumpt = 30:30:360;
% if the keyword is differential pressure set it to ISIP. For dP and
% FlowRate adjust the nparam value. Check that the parameters are
% valid
nparam_orig = nparam;
keyword_orig = keyword;
switch upper(keyword)
case 'DP'
keyword = 'ISIP';
if isempty(find(dP == nparam, 1))
error('Parameter value invalid');
end
nparam = (nparam + 36.824) * 1e6;
case 'FLOWRATE'
if isempty(find(flow == nparam, 1))
error('Parameter value invalid');
end
nparam = round((nparam * 86400)/100)*100;
case 'PUMPTIME'
if isempty(find(pumpt == nparam, 1))
error('Parameter value invalid');
end
end
% generate empty matrices for the lateral distance calculation
all_x = [];
all_y = [];
all_z = [];
all_stress = [];
% have to get the data for high and low res for the distance
% calculation, will only plot the required one.
for interval = [20 50]
%set the index
if interval == 20
ind = 1;
else
ind = 2;
end
% generate the filenames of each of the mechanisms for 100% i, 100% s
% and 100% d
f_i = fullfile(root_folder, ...
[keyword '_' num2str(nparam) '_i100-s0-d0_' ...
num2str(interval) 'm.mat']);
f_s = fullfile(root_folder, ...
[keyword '_' num2str(nparam) '_i0-s100-d0_' ...
num2str(interval) 'm.mat']);
f_d = fullfile(root_folder, ...
[keyword '_' num2str(nparam) '_i0-s0-d100_' ...
num2str(interval) 'm.mat']);
% check the mat files are there
if ~exist(f_i, 'file') || ~exist(f_d, 'file') || ~exist(f_s, 'file')
error(['Not all mat files are present: ' num2str(nparam_orig)]);
end
% Get the data from the mat files
m_i = load(f_i);
m_s = load(f_s);
m_d = load(f_d);
% save the x, y, z
x_i = m_i.x;
y_i = m_i.y;
z_i = m_i.z;
% x = x_i(:);
% y = y_i(:);
zz = zeros(length(x_i), 50);
% get the correct stress matrix
switch (upper(stresstype))
case 'C'
c_i = weightings(1) .* m_i.coulomb;
c_s = weightings(2) .* m_s.coulomb;
c_d = weightings(3) .* m_d.coulomb;
st = 'Coulomb';
case 'N'
c_i = weightings(1) .* m_i.normal;
c_s = weightings(2) .* m_s.normal;
c_d = weightings(3) .* m_d.normal;
st = 'Normal';
case 'S'
c_i = weightings(1) .* m_i.shear;
c_s = weightings(2) .* m_s.shear;
c_d = weightings(3) .* m_d.shear;
st = 'Shear';
otherwise
error('Stress type not recognised');
end
% Add the stresses together
c = c_i + c_s + c_d;
%preallocate the maximum stress matrix
lux = length(unique(x_i));
luy = length(unique(y_i));
stress_m = zeros(lux, luy, 50);
s_m = zeros(lux*luy, 50);
%xyz = zeros(50,1);
% get the max stress for the required monte carlo iteration
for jj = 1:50
[stress_m(:,:,jj), xyz(jj)] = max_stress(x_i, y_i, z_i, c(:,:,jj)); %#ok<AGROW>
s_m(:,jj) = reshape(stress_m(:,:,jj), [], 1);
zz(:,jj) = xyz(jj).z(:);
end
all_x = [all_x; x_i(:)];
all_y = [all_y; y_i(:)];
all_z = [all_z; zz];
all_stress = [all_stress; s_m];
% plot the stress map
if resolution == interval% get the fractures
fx = m_i.fracturesX;
fy = m_i.fracturesY;
fracX = fx{mc};
fracY = fy{mc};
fig = plot_stress(fracX, fracY, stress_m(:,:,mc), xyz(mc), plot_fracs, st, resolution);
end
end
% calculate the lateral distance
[mean_ld, median_d, max_ld, min_ld] = ...
lateralDistance(all_x, all_y, all_stress, ldvalue);
% add the text with poisson ratio, youngs modulus, regional stress,
% flow, area and lateral distance
filename = [keyword_orig '_' st '_' num2str(resolution) 'm_' ...
num2str(nparam_orig) '_' ...
num2str(weightings(1)*100) 'i' num2str(weightings(2)*100) 's' num2str(weightings(3)*100) 'r__' ...
num2str(mc)];
saveText(filename, m_i, keyword, [mean_ld, median_d, max_ld, min_ld], ldvalue, mc);
saveas(fig, [filename '.jpg']);
end
function saveText(filename, m_i, keyword, latdist, latdistthreshold, mc)
% Function which generates additional information
% and returns it in a cell array.
fid = fopen([filename '.txt'], 'w');
pr = m_i.poisson;
ym = m_i.youngs * 100000 / 1e9;
rs = m_i.regional_stress;
fric = m_i.friction;
dimensions = m_i.dimensions;
res = dimensions(3);
param_id = m_i.param_id;
switch keyword
case 'ISIP'
keyword = 'dP';
param = (param_id * 1e-6) - 36.824;
param2 = '';
unit = 'MPa';
case 'FlowRate'
param = param_id;
unit = 'm^3/d';
unit2 = 'm^3/s';
param2 = [' (' num2str(param/86400, 5) ' ' unit2 ')'];
case 'PumpTime'
param = param_id;
unit = 'mins';
param2 = '';
end
fprintf(fid, 'Calculation for: \n');
fprintf(fid, '%s %u %s %s \n%u m resolution \n', keyword, param, unit, param2);
fprintf(fid, 'Monte Carlo Iteration %u\n \n', mc);
fprintf(fid, 'Geological Conditions: \n');
fprintf(fid, 'Poisson ratio %2.2f\n', pr);
fprintf(fid, 'Young''s modulus %2.2f GPa\n', ym);
fprintf(fid, 'Friction %2.2f\n\n', fric);
fprintf(fid, 'Regional Stress: %s \n\n', mat2str(rs));
fprintf(fid, 'Lateral Respect Distance (lrd) data %2.3f MPa threshold\n', latdistthreshold);
fprintf(fid, 'Mean lrd %4.2f \n', latdist(1)*1000);
fprintf(fid, 'Median lrd %4.2f \n', latdist(2)*1000);
fprintf(fid, 'Maximum lrd %4.2f \n', latdist(3)*1000);
fprintf(fid, 'Minimum lrd %4.2f \n', latdist(4)*1000);
fclose(fid);
end
function [stress_max, xyz] = max_stress(x, y, z, stress)
% MAX_STRESS compute the max or min stress over all depths
% [MS, XYZ] = MAX_STRESS(X, Y, Z, STRESS)
% Takes the X, Y, Z and STRESS vectors, representing a grid of
% points and stress and returns the maximum and minimum stresses at
% each radius in MPa
stress = stress/10; % divide by 10 to get MPa
% reshape the data to be a cube
stress = reshape(stress, length(unique(x)), length(unique(y)), ...
length(unique(z)), []);
stress_max = stress;
xyz.x = x(:);
xyz.y = y(:);
xyz.z = z(:);
if ndims(stress_max) > 2 %#ok<ISMAT> % check how many levels of z we have
% calculate using abs max at each x,y and make sign right
[stress_max, max_position] = max(abs(stress),[],3);
[~, min_position] = min(stress,[],3);
stress_sign = (min_position == max_position) .* -1.0 ;
stress_sign(stress_sign == 0) = 1;
stress_max = stress_max .* stress_sign;
xyz.z = z(max_position);
end
end
function fig = plot_stress(fracX, fracY, stress, xyz, plot_fracs, stresstype, resolution)
%PLOT_STRESS. Plots the stress maps
% PLOT_STRESS(FRACX, FRACY, STRESS, XYZ, PLOT_FRACS)
% Reads in the maximum STRESS and XYZ and plots the stress maps, with
% the option to overlay the fractures
%
% FRACX and FRACY are the x and y coordinates of the fractures
% respectively.
% STRESS is the stress in a m x n matrix
% XYZ is a structure containing the x, y and z values as vectors
% PLOT_FRACS is logical and states whether the fractures should be
% plotted or not.
if ~exist('plot_fracs', 'var') || isempty(plot_fracs)
plot_fracs = 1;
end
% set the logical for plotting the fractures and the fontsize
PLOT_FRACS = logical(plot_fracs);
FONTSIZE = 14;
% set the boundaries for clipping the data
if resolution == 50
clip = [-4e-4 4e-4];
interval = 1e-5;
else
clip = [-4e-5 4e-5];
interval = 1e-6;
end
% clip the data
stress(stress < clip(1)) = clip(1);
stress(stress > clip(2)) = clip(2);
% set the colour levels
c_levels = clip(1):interval:clip(2);
% generate a new figure
scrsz = get(groot,'ScreenSize');
fig = figure('Position',[10 scrsz(4)/2-100 scrsz(3)/2 scrsz(4)/2]);
% set the colourmap
colormap(anatolia);
% contour plot the stress and set the fontsize of the axis
contourf(unique(xyz.x), unique(xyz.y), stress, c_levels, 'LineColor', 'none');
shading interp
set(gca, 'FontSize', FONTSIZE);%ax.FontSize = FONTSIZE;
% plot contour lines - black for positive and grey for negative
hold on
contour(unique(xyz.x), unique(xyz.y), stress, 0:1e-5:1e-4, 'k');
contour(unique(xyz.x), unique(xyz.y), stress, 0:-1e-5:-1e-4, 'LineColor', [.5 .5 .5]);
hold off
% set the labels
xlabel('X (m)','FontSize',FONTSIZE)
ylabel('Y (m)','FontSize',FONTSIZE)
caxis (clip)
% set axis details
axis equal
axis([min(xyz.x) max(xyz.x) min(xyz.y) max(xyz.y)]);
if PLOT_FRACS
fcolour = [0 .7 0]; % green
% get the number of fractures
m = length(fracX(1,:));
% cycle round and plot
for k = 1:m
line(fracX(:,k), fracY(:,k), 'LineStyle', '-', 'Color', fcolour);
end
end
% add the colourbar
if resolution == 50
hcb = colorbar('Location', 'EastOutside',...
'Ticks', [-4e-4,-2e-4,0,2e-4,4e-4],...
'TickLabels', {'-0.0004', '-0.0002', '0' '0.0002', '0.0004'}, ...
'FontSize', 14);
else
hcb = colorbar('Location', 'EastOutside',...
'Ticks', [-4e-5,-2e-5,0,2e-5,4e-5],...
'TickLabels', {'-0.00004', '-0.00002', '0' '0.00002', '0.00004'}, ...
'FontSize', 14);
end
title(hcb, [stresstype ' Stress (MPa)']);
%set(get(get(hcb, 'Label'), 'String'), [stresstype ' Stress (MPa)']);
% set the renderer
set(gcf,'renderer','painters')
end
function [a,b,r,s] = conv_hull(stress, x, y)
%CONV_HULL. Obtains all the points of the max stress v radius with the conve hull
%bound
%
% [A, B, R, S] = CONV_HULL(STRESS, X, Y)
%
% STRESS is the stress in a mxn matrix
% X is the x values as vectors
% Y is the y values as vectors
%
% A and B are the parameters to define the convex hull
% R is the radius in km and S is the max stress
% set the centre point
pt = [0,0];
% calculate the radius from the centre of the fractures and sort
radius = sqrt((x-pt(1)).^2 + (y-pt(2)).^2);
[radius, idx]= sort(radius);
% sort the stress into the same order as the radius
stress = stress(:);
stress = stress(idx);
% get the unique radius
r = unique(radius);
% create a new array for the max stress at each radius
s = zeros(size(r));
% loop round to populate s
for n = 1:length(s)
ind = radius == r(n);
s(n) = max(stress(ind));
end
% find all values of stress less than 0 and get rid of them
ind = s <= 0;
r(ind) = [];
s(ind) = [];
% Get rid of any radius close in (<0.25 km)
r1 = r(r>0.25);
s1 = s(r>0.25);
% generate the convexhull
ind = convhull(log(r1), log(s1));
% number of points
i = numel(ind);
% loop round to find the first point on the curve
while (r1(ind(i))<0.750) %& i>1
i = i-1;
end
R1 = r1(ind(i));
E1 = s1(ind(i));
R2 = r1(ind(i+1));
E2 = s1(ind(i+1));
b = log(E1/E2)/log(R1/R2);
a = E1/R1^b;
end
function descr = outputText(m_i, keyword, latdist, latdistthreshold, mc)
% Function which generates additional information
% and returns it in a cell array.
pr = m_i.poisson;
ym = m_i.youngs * 100000 / 1e9;
rs = m_i.regional_stress;
fric = m_i.friction;
dimensions = m_i.dimensions;
res = dimensions(3);
param_id = m_i.param_id;
switch keyword
case 'ISIP'
keyword = 'dP';
param = (param_id * 1e-6) - 36.824;
param2 = '';
unit = 'MPa';
case 'FlowRate'
param = param_id;
unit = 'm^3/d';
unit2 = 'm^3/s';
param2 = [' (' num2str(param/86400, 5) ' ' unit2 ')'];
case 'PumpTime'
param = param_id;
unit = 'mins';
param2 = '';
end
descr = {'Calculation for:';
[keyword ': ' num2str(param, 5) ' ' unit param2 ', ' ...
num2str(res) 'm resolution'];
['Monte Carlo iteration: ' num2str(mc)];
' ';
'Geological Conditions:';
['Poisson ratio ' num2str(pr)];
['Young''s modulus ' num2str(ym) ' GPa'];
['Friction ' num2str(fric)]
'';
'Regional Stress:';
mat2str(rs)
'';
['Lateral Respect Distance (lrd) data ' num2str(latdistthreshold) ...
' MPa threshold: ']
['Mean lrd ' num2str(latdist(1)*1000) ' m']
['Median lrd ' num2str(latdist(2)*1000) ' m']
['Maximum lrd ' num2str(latdist(3)*1000) ' m']
['Minimum lrd ' num2str(latdist(1)*1000) ' m'];
};
end
function [mean_ld, median_ld, max_ld, min_ld] = lateralDistance(x, y, stress, latdistthreshold)
%LATERALDISTANCE Calculates the lateral distance using a convex
% hull
%
% [MEAN_LD, MEDIAN_LD, MAX_LD, MIN_LD] = LATERALDISTANCE(X, Y, Z, STRESS, FRAC_CENT)
%
% X is the x values as a vector
% Y is the y values as a vector
% STRESS is the stress in a m x n matrix
% LATDISTTHRESHOLD is the stress threshold to calculate the lateral
% respect distance for
% check if we are in km, if not, convert
if max(abs(x)) > 10
x = x/1000;
y = y/1000;
end
% create a radius vector
radius = zeros(1,50);
% loop round all monte carlo iterations
for jj = 1:50
[a, b] = conv_hull(stress(:,jj), x, y);
% compute the radius
radius(jj) = (latdistthreshold/a)^(1/b);
end
mean_ld = mean(radius);
median_ld = median(radius);
max_ld = max(radius);
min_ld = min(radius);
end
function cmap = anatolia()
% the Anatolia colormap
cmap = ...
[0.21176470816135406 0.15294118225574493 0.99215686321258545;
0.21398124098777771 0.1877237856388092 0.99232739210128784;
0.21619778871536255 0.22250640392303467 0.99249786138534546;
0.21841432154178619 0.25728902220726013 0.99266839027404785;
0.22063086926937103 0.29207161068916321 0.99283885955810547;
0.22284740209579468 0.32685422897338867 0.99300938844680786;
0.22506394982337952 0.36163684725761414 0.99317985773086548;
0.22728048264980316 0.39641943573951721 0.99335038661956787;
0.22949701547622681 0.43120205402374268 0.99352091550827026;
0.23171356320381165 0.46598467230796814 0.99369138479232788;
0.23393009603023529 0.5007672905921936 0.99386191368103027;
0.23614664375782013 0.53554987907409668 0.99403238296508789;
0.23836317658424377 0.57033246755599976 0.99420291185379028;
0.24057972431182861 0.60511511564254761 0.9943733811378479;
0.24279625713825226 0.63989770412445068 0.99454391002655029;
0.2450128048658371 0.67468029260635376 0.99471437931060791;
0.24722933769226074 0.70946294069290161 0.9948849081993103;
0.24944587051868439 0.74424552917480469 0.9950554370880127;
0.25166240334510803 0.77902811765670776 0.99522590637207031;
0.25387895107269287 0.81381076574325562 0.99539643526077271;
0.25609549880027771 0.84859335422515869 0.99556690454483032;
0.25831204652786255 0.88337594270706177 0.99573743343353271;
0.260528564453125 0.91815859079360962 0.99590790271759033;
0.26274511218070984 0.9529411792755127 0.99607843160629272;
0.35490196943283081 0.958823561668396 0.99656862020492554;
0.44705882668495178 0.96470588445663452 0.99705880880355835;
0.53921568393707275 0.970588207244873 0.99754899740219116;
0.63137257099151611 0.97647058963775635 0.99803924560546875;
0.7235293984413147 0.98235297203063965 0.99852943420410156;
0.81568628549575806 0.98823529481887817 0.99901962280273438;
0.90784311294555664 0.9941176176071167 0.99950981140136719;
1 1 1;
1 1 1;
0.99901962280273438 0.99215686321258545 0.91372549533843994;
0.99803924560546875 0.9843137264251709 0.82745099067687988;
0.99705880880355835 0.97647058963775635 0.74117648601531982;
0.99607843160629272 0.9686274528503418 0.65490198135375977;
0.9950980544090271 0.96078431606292725 0.56862747669219971;
0.9941176176071167 0.9529411792755127 0.48235294222831726;
0.99313724040985107 0.94509804248809814 0.3960784375667572;
0.99215686321258545 0.93725490570068359 0.30980393290519714;
0.989599347114563 0.89650470018386841 0.29991474747657776;
0.98704177141189575 0.85575449466705322 0.29002559185028076;
0.98448425531387329 0.815004289150238 0.28013640642166138;
0.981926679611206 0.77425402402877808 0.27024725079536438;
0.97936916351318359 0.73350381851196289 0.260358065366745;
0.97681158781051636 0.69275361299514771 0.25046887993812561;
0.9742540717124939 0.65200340747833252 0.24057972431182861;
0.97169649600982666 0.61125320196151733 0.23069053888320923;
0.9691389799118042 0.57050299644470215 0.22080136835575104;
0.966581404209137 0.529752790927887 0.21091219782829285;
0.9640238881111145 0.48900255560874939 0.20102302730083466;
0.96146631240844727 0.4482523500919342 0.19113385677337646;
0.9589087963104248 0.407502144575119 0.18124467134475708;
0.95635122060775757 0.36675190925598145 0.17135550081729889;
0.95379370450973511 0.32600170373916626 0.1614663302898407;
0.95123612880706787 0.28525149822235107 0.15157715976238251;
0.94867861270904541 0.24450127780437469 0.14168798923492432;
0.94612103700637817 0.20375107228755951 0.13179880380630493;
0.94356352090835571 0.16300085186958313 0.12190964072942734;
0.94100594520568848 0.12225063890218735 0.11202046275138855;
0.938448429107666 0.081500425934791565 0.10213129222393036;
0.93589085340499878 0.040750212967395782 0.092242114245891571;
0.93333333730697632 0 0.08235294371843338];
end