% Script for Wellbore Bridging Simulation % Author: Oluwafemi Oyedokun % Reservoir Fluid Type: Oil and Gas % Objective: Wellbore Bridging Only % Assumptions % 1.Plane Strain approimation is assumed. % 2.Infinite acting reservoir(s) (i.e transient flow) % 3. Tranversely isotropic producing formations (affecting permeability % only) % 4. No skin factor due to formation damage, except well-reservoir geometry clc clear % Loading the Comsol Model Object BoraSim into MATLAB import com.comsol.model.* import com.comsol.model.util.* model=mphload('BoraSim-I.mph'); %model=ModelUtil.model('Model'); % model.study('std1').run; model.sol('sol2').run; %% %DATA LINKED TO COMSOL DESKTOP %Importing the wellbore geometry data from Comsol Desktop to Matlab aa = mphevaluate(model,'dcov')*ones(15,1); %The first fifteen openhole diameters in Vertical Section in Comsol base unit(meters) ab = mphevaluate(model,'dcoc')*ones(3,1);% The diameters of next three openhole sections in curved section in comsol base unit(meters) ac = mphevaluate(model,'dcol'); %Diameters of Lateral Openhole Sections in comsol base unit(meters) a = 0.5*39.3701*[aa(1) aa(2) aa(3) aa(4) aa(5) aa(6) aa(7) aa(8) aa(9) aa(10)... aa(11) aa(12) aa(13) aa(14) aa(15) ab(1) ab(2) ab(3) ac ac ac... ac ac ac ac ac ac ac]; %Initial openhole radius before breakout in inches ba = 0.5*mphevaluate(model,'dr'); %Radius of the riser in Comsol base unit(meters) bb = 0.5*mphevaluate(model,'dcv');%Inner radius of Last casing in the vertical section in comsol unit(meters) bc = 0.5*mphevaluate(model,'dLv');%Inner radius of Liner in the vertical section in comsol unit(meters) bd = 0.5*mphevaluate(model,'dcc');%Inner radius of casing in the curved section in comsol unit(meters) be = 0.5*mphevaluate(model,'dcl');%Inner radius of casing in the lateral section in comsol unit(meters) rc = 39.3701*[ba bb bc bd be]; %Internal radii of cased and liner sections respectively, in inches ca = mphevaluate(model,'hw');%water depth or riser length in comsol unit(meters) cb = mphevaluate(model,'csv');%casing length in the vertical section in comsol unit(meters) cc = mphevaluate(model,'hLv');%liner length in the vertical section in comsol unit(meters) Inc = pi/180*mphevaluate(model,'Incc');%Curved-cased section run angle in rad Rcc = mphevaluate(model,'Rc');%Radius of curvature of the curved section in meters cd = pi/180*Rcc*Inc;% Casing run length in the curved section in meters ce = mphevaluate(model,'Lc');%Casing length in the lateral section in meters cf = mphevaluate(model,'ucsv');%Openhole Length in the vertical section of the wellbore hc = 3.28084*[ca cb cc cd ce];%Differential Length of casing and liner sections in feet %Importing the Width of Artificial Wellbore Enlargement in meters from COMSOL tL1 = mphevaluate(model,'tL1');tL2 = mphevaluate(model,'tL2');tL3 = mphevaluate(model,'tL3'); tL4 = mphevaluate(model,'tL4');tL5 = mphevaluate(model,'tL5');tL6 = mphevaluate(model,'tL6'); tL7 = mphevaluate(model,'tL7');tL8 = mphevaluate(model,'tL8');tL9 = mphevaluate(model,'tL9'); tL10 = mphevaluate(model,'tL10');tL11 = mphevaluate(model,'tL11');tL12 = mphevaluate(model,'tL12'); tL13 = mphevaluate(model,'tL13');tL14 = mphevaluate(model,'tL14');tL15 = mphevaluate(model,'tL15'); tL16 = mphevaluate(model,'tL16');tL17 = mphevaluate(model,'tL17');tL18 = mphevaluate(model,'tL18'); tL19 = mphevaluate(model,'tL19');tL20 = mphevaluate(model,'tL20');tL21 = mphevaluate(model,'tL21'); tL22 = mphevaluate(model,'tL22');tL23 = mphevaluate(model,'tL23');tL24 = mphevaluate(model,'tL24'); tL25 = mphevaluate(model,'tL25');tL26 = mphevaluate(model,'tL26');tL27 = mphevaluate(model,'tL27'); tL28 = mphevaluate(model,'tL28'); rbkt_int = 39.3701*[tL1 tL2 tL3 tL4 tL5 tL6 tL7 tL8 tL9 tL10 tL11 tL12 tL13 tL14... tL15 tL16 tL17 tL18 tL19 tL20 tL21 tL22 tL23 tL24 tL25 tL26 tL27 tL28]; % Artificial Breakout in inches %Openhole Well Length in Comsol Base Unit(meters) h1 = mphevaluate(model,'h1v'); h2 = mphevaluate(model,'h2v'); h3 = mphevaluate(model,'h3v'); h4 = mphevaluate(model,'h4v'); h5 = mphevaluate(model,'h5v'); h6 = mphevaluate(model,'h6v'); h7 = mphevaluate(model,'h7v'); h8 = mphevaluate(model,'h8v'); h9 = mphevaluate(model,'h1v'); h10 = mphevaluate(model,'h10v'); h11 = mphevaluate(model,'h11v'); h12 = mphevaluate(model,'h12v'); h13 = mphevaluate(model,'h13v'); h14 = mphevaluate(model,'h14v'); h15 = mphevaluate(model,'h15v'); h16 = mphevaluate(model,'h1c'); h17 = mphevaluate(model,'h2c'); h18 = mphevaluate(model,'h3c');h19 = mphevaluate(model,'hL1'); h20 = mphevaluate(model,'hL2'); h21 = mphevaluate(model,'hL3'); h22 = mphevaluate(model,'hL4'); h23 = mphevaluate(model,'hL5'); h24 = mphevaluate(model,'hL6'); h25 = mphevaluate(model,'hL7'); h26 = mphevaluate(model,'hL8'); h27 = mphevaluate(model,'hL9'); h28 = mphevaluate(model,'hL10'); h=3.28084*[h1 h2 h3 h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16 h17 h18 h19 h20... h21 h22 h23 h24 h25 h26 h27 h28]; %Openhole Section Length in feet %Well Angles from Comsol in Rad Inc1 = mphevaluate(model,'Inc1c'); Inc2 = mphevaluate(model,'Inc2c'); Inc3 = mphevaluate(model,'Inco'); %Other Parameters alpL = mphevaluate(model,'alpLat'); kLat = mphevaluate(model,'kLat'); % Fluid Flow Parameters Extracted at Mid Points of Each Openhole Layer %Extracting the velocity and pressure values at the following coordinates: %e.g x12= 1-Layer 1 2-Point 2 of Layer 1; y1-Mid depth of Layer 1 from origin %Midline Coordinates in the vertical section of the wellbore x11=-aa(1)/4;x12=0;x13=aa(1)/4;x14=-(2/3*tL1+aa(1)/2);x15=-(1/3*tL1+aa(1)/2);x16=(1/3*tL1+aa(1)/2);x17=(2/3*tL1+aa(1)/2); y1=-(ca/2+cb+cc+h1/2); x21=-aa(2)/4;x22=0;x23=aa(2)/4;x24=-(2/3*tL2+aa(2)/2);x25=-(1/3*tL2+aa(2)/2);x26=(1/3*tL2+aa(2)/2);x27=(2/3*tL2+aa(2)/2); y2=-(ca/2+cb+cc+h1+h2/2); x31=-aa(3)/4;x32=0;x33=aa(3)/4;x34=-(2/3*tL3+aa(3)/2);x35=-(1/3*tL3+aa(3)/2);x36=(1/3*tL3+aa(3)/2);x37=(2/3*tL3+aa(3)/2); y3=-(ca/2+cb+cc+h1+h2+h3/2); x41=-aa(4)/4;x42=0;x43=aa(4)/4;x44=-(2/3*tL4+aa(4)/2);x45=-(1/3*tL4+aa(4)/2);x46=(1/3*tL4+aa(4)/2);x47=(2/3*tL4+aa(4)/2); y4=-(ca/2+cb+cc+h1+h2+h3+h4/2); x51=-aa(5)/4;x52=0;x53=aa(5)/4;x54=-(2/3*tL5+aa(5)/2);x55=-(1/3*tL5+aa(5)/2);x56=(1/3*tL5+aa(5)/2);x57=(2/3*tL5+aa(5)/2); y5=-(ca/2+cb+cc+h1+h2+h3+h4+h5/2); x61=-aa(6)/4;x62=0;x63=aa(6)/4;x64=-(2/3*tL6+aa(6)/2);x65=-(1/3*tL6+aa(6)/2);x66=(1/3*tL6+aa(6)/2);x67=(2/3*tL6+aa(6)/2); y6=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6/2); x71=-aa(7)/4;x72=0;x73=aa(7)/4;x74=-(2/3*tL7+aa(7)/2);x75=-(1/3*tL7+aa(7)/2);x76=(1/3*tL7+aa(7)/2);x77=(2/3*tL7+aa(7)/2); y7=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7/2); x81=-aa(8)/4;x82=0;x83=aa(8)/4;x84=-(2/3*tL8+aa(8)/2);x85=-(1/3*tL8+aa(8)/2);x86=(1/3*tL8+aa(8)/2);x87=(2/3*tL8+aa(8)/2); y8=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8/2); x91=-aa(9)/4;x92=0;x93=aa(9)/4;x94=-(2/3*tL9+aa(9)/2);x95=-(1/3*tL9+aa(9)/2);x96=(1/3*tL9+aa(9)/2);x97=(2/3*tL9+aa(9)/2); y9=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9/2); x101=-aa(10)/4;x102=0;x103=aa(10)/4;x104=-(2/3*tL10+aa(10)/2);x105=-(1/3*tL10+aa(10)/2);x106=(1/3*tL10+aa(10)/2);x107=(2/3*tL10+aa(10)/2); y10=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9+h10/2); x111=-aa(11)/4;x112=0;x113=aa(11)/4;x114=-(2/3*tL11+aa(11)/2);x115=-(1/3*tL11+aa(11)/2);x116=(1/3*tL11+aa(11)/2);x117=(2/3*tL11+aa(11)/2); y11=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9+h10+h11/2); x121=-aa(12)/4;x122=0;x123=aa(12)/4;x124=-(2/3*tL12+aa(12)/2);x125=-(1/3*tL12+aa(12)/2);x126=(1/3*tL12+aa(12)/2);x127=(2/3*tL12+aa(12)/2); y12=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9+h10+h11+h12/2); x131=-aa(13)/4;x132=0;x133=aa(13)/4;x134=-(2/3*tL13+aa(13)/2);x135=-(1/3*tL13+aa(13)/2);x136=(1/3*tL13+aa(13)/2);x137=(2/3*tL13+aa(13)/2); y13=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9+h10+h11+h12+h13/2); x141=-aa(14)/4;x142=0;x143=aa(14)/4;x144=-(2/3*tL14+aa(14)/2);x145=-(1/3*tL14+aa(14)/2);x146=(1/3*tL14+aa(14)/2);x147=(2/3*tL14+aa(14)/2); y14=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9+h10+h11+h12+h13+h14/2); x151=-aa(15)/4;x152=0;x153=aa(15)/4;x154=-(2/3*tL15+aa(15)/2);x155=-(1/3*tL15+aa(15)/2);x156=(1/3*tL15+aa(15)/2);x157=(2/3*tL15+aa(15)/2); y15=-(ca/2+cb+cc+h1+h2+h3+h4+h5+h6+h7+h8+h9+h10+h11+h12+h13+h14+h15/2); %Midline Coordinates in the Curved Section of the Wellbore Inchlf1 = 1/2*(Inc+Inc1); Inchlf2 = 1/2*(Inc1+Inc2); Inchlf3 = 1/2*(Inc2+Inc3); x161=-ab(1)/4+(Rcc+ab(1)/4)*(1-cos(Inchlf1));y161=-(ca/2+cb+cc+cf)-(Rcc+ab(1)/4)*sin(Inchlf1); x162=Rcc*(1-cos(Inchlf1));y162=-(ca/2+cb+cc+cf)-Rcc*sin(Inchlf1); x163=ab(1)/4+(Rcc-ab(1)/4)*(1-cos(Inchlf1));y163=-(ca/2+cb+cc+cf)-(Rcc-ab(1)/4)*sin(Inchlf1); x164=-(2/3*tL16+ab(1)/2)+(Rcc+2/3*tL16+ab(1)/2)*(1-cos(Inchlf1));y164=-(ca/2+cb+cc+cf)-(Rcc+2/3*tL16+ab(1)/2)*sin(Inchlf1); x165=-(1/3*tL16+ab(1)/2)+(Rcc+1/3*tL16+ab(1)/2)*(1-cos(Inchlf1));y165=-(ca/2+cb+cc+cf)-(Rcc+1/3*tL16+ab(1)/2)*sin(Inchlf1); x166=(1/3*tL16+ab(1)/2)+(Rcc-1/3*tL16-ab(1)/2)*(1-cos(Inchlf1));y166=-(ca/2+cb+cc+cf)-(Rcc-1/3*tL16-ab(1)/2)*sin(Inchlf1); x167=(2/3*tL16+ab(1)/2)+(Rcc-2/3*tL16-ab(1)/2)*(1-cos(Inchlf1));y167=-(ca/2+cb+cc+cf)-(Rcc-2/3*tL16-ab(1)/2)*sin(Inchlf1); x171=-ab(2)/4+(Rcc+ab(2)/4)*(1-cos(Inchlf2));y171=-(ca/2+cb+cc+cf)-(Rcc+ab(2)/4)*sin(Inchlf2); x172=Rcc*(1-cos(Inchlf2));y172=-(ca/2+cb+cc+cf)-Rcc*sin(Inchlf2); x173=ab(2)/4+(Rcc-ab(2)/4)*(1-cos(Inchlf2));y173=-(ca/2+cb+cc+cf)-(Rcc-ab(2)/4)*sin(Inchlf2); x174=-(2/3*tL17+ab(2)/2)+(Rcc+2/3*tL17+ab(2)/2)*(1-cos(Inchlf2));y174=-(ca/2+cb+cc+cf)-(Rcc+2/3*tL17+ab(2)/2)*sin(Inchlf2); x175=-(1/3*tL17+ab(2)/2)+(Rcc+1/3*tL17+ab(2)/2)*(1-cos(Inchlf2));y175=-(ca/2+cb+cc+cf)-(Rcc+1/3*tL17+ab(2)/2)*sin(Inchlf2); x176=(1/3*tL17+ab(2)/2)+(Rcc-1/3*tL17-ab(2)/2)*(1-cos(Inchlf2));y176=-(ca/2+cb+cc+cf)-(Rcc-1/3*tL17-ab(2)/2)*sin(Inchlf2); x177=(2/3*tL17+ab(2)/2)+(Rcc-2/3*tL17-ab(2)/2)*(1-cos(Inchlf2));y177=-(ca/2+cb+cc+cf)-(Rcc-2/3*tL17-ab(2)/2)*sin(Inchlf2); x181=-ab(3)/4+(Rcc+ab(3)/4)*(1-cos(Inchlf3));y181=-(ca/2+cb+cc+cf)-(Rcc+ab(3)/4)*sin(Inchlf3); x182=Rcc*(1-cos(Inchlf3));y182=-(ca/2+cb+cc+cf)-Rcc*sin(Inchlf3); x183=ab(3)/4+(Rcc-ab(2)/4)*(1-cos(Inchlf3));y183=-(ca/2+cb+cc+cf)-(Rcc-ab(3)/4)*sin(Inchlf3); x184=-(2/3*tL18+ab(3)/2)+(Rcc+2/3*tL18+ab(3)/2)*(1-cos(Inchlf3));y184=-(ca/2+cb+cc+cf)-(Rcc+2/3*tL18+ab(3)/2)*sin(Inchlf3); x185=-(1/3*tL18+ab(3)/2)+(Rcc+1/3*tL18+ab(3)/2)*(1-cos(Inchlf3));y185=-(ca/2+cb+cc+cf)-(Rcc+1/3*tL18+ab(3)/2)*sin(Inchlf3); x186=(1/3*tL18+ab(3)/2)+(Rcc-1/3*tL18-ab(3)/2)*(1-cos(Inchlf3));y186=-(ca/2+cb+cc+cf)-(Rcc-1/3*tL18-ab(3)/2)*sin(Inchlf3); x187=(2/3*tL18+ab(3)/2)+(Rcc-2/3*tL18-ab(3)/2)*(1-cos(Inchlf3));y187=-(ca/2+cb+cc+cf)-(Rcc-2/3*tL18-ab(3)/2)*sin(Inchlf3); %Midline Coordinates in the Lateral Section of the Wellbore x192 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180); y192 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180); x191 = x192-ac/4*cos(pi/180*Inc3);y191 = y192-ac/4*sin(pi/180*Inc3); x193 = x192+ac/4*cos(pi/180*Inc3);y193 = y192+ac/4*sin(pi/180*Inc3); x194 = x192-(ac/2+2/3*tL19)*cos(pi/180*Inc3);y194 = y192-(ac/2+2/3*tL19)*sin(pi/180*Inc3); x195 = x192-(ac/2+1/3*tL19)*cos(pi/180*Inc3);y195 = y192-(ac/2+1/3*tL19)*sin(pi/180*Inc3); x196 = x192+(ac/2+1/3*tL19)*cos(pi/180*Inc3);y196 = y192+(ac/2+1/3*tL19)*sin(pi/180*Inc3); x197 = x192+(ac/2+2/3*tL19)*cos(pi/180*Inc3);y197 = y192+(ac/2+2/3*tL19)*sin(pi/180*Inc3); x202 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+h20)*sin(Inc3*pi/180); y202 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+h20)*cos(Inc3*pi/180); x201 = x202-ac/4*cos(pi/180*Inc3);y201 = y202-ac/4*sin(pi/180*Inc3); x203 = x202+ac/4*cos(pi/180*Inc3);y203 = y202+ac/4*sin(pi/180*Inc3); x204 = x202-(ac/2+2/3*tL20)*cos(pi/180*Inc3);y204 = y202-(ac/2+2/3*tL20)*sin(pi/180*Inc3); x205 = x202-(ac/2+1/3*tL20)*cos(pi/180*Inc3);y205 = y202-(ac/2+1/3*tL20)*sin(pi/180*Inc3); x206 = x202+(ac/2+1/3*tL20)*cos(pi/180*Inc3);y206 = y202+(ac/2+1/3*tL20)*sin(pi/180*Inc3); x207 = x202+(ac/2+2/3*tL20)*cos(pi/180*Inc3);y207 = y202+(ac/2+2/3*tL20)*sin(pi/180*Inc3); x212 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+h21)*sin(Inc3*pi/180); y212 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+h21)*cos(Inc3*pi/180); x211 = x212-ac/4*cos(pi/180*Inc3);y211 = y212-ac/4*sin(pi/180*Inc3); x213 = x212+ac/4*cos(pi/180*Inc3);y213 = y212+ac/4*sin(pi/180*Inc3); x214 = x212-(ac/2+2/3*tL21)*cos(pi/180*Inc3);y214 = y212-(ac/2+2/3*tL21)*sin(pi/180*Inc3); x215 = x212-(ac/2+1/3*tL21)*cos(pi/180*Inc3);y215 = y212-(ac/2+1/3*tL21)*sin(pi/180*Inc3); x216 = x212+(ac/2+1/3*tL21)*cos(pi/180*Inc3);y216 = y212+(ac/2+1/3*tL21)*sin(pi/180*Inc3); x217 = x212+(ac/2+2/3*tL21)*cos(pi/180*Inc3);y217 = y212+(ac/2+2/3*tL21)*sin(pi/180*Inc3); x222 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+h22)*sin(Inc3*pi/180); y222 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+h22)*cos(Inc3*pi/180); x221 = x222-ac/4*cos(pi/180*Inc3);y221 = y222-ac/4*sin(pi/180*Inc3); x223 = x222+ac/4*cos(pi/180*Inc3);y223 = y222+ac/4*sin(pi/180*Inc3); x224 = x222-(ac/2+2/3*tL22)*cos(pi/180*Inc3);y224 = y222-(ac/2+2/3*tL22)*sin(pi/180*Inc3); x225 = x222-(ac/2+1/3*tL22)*cos(pi/180*Inc3);y225 = y222-(ac/2+1/3*tL22)*sin(pi/180*Inc3); x226 = x222+(ac/2+1/3*tL22)*cos(pi/180*Inc3);y226 = y222+(ac/2+1/3*tL22)*sin(pi/180*Inc3); x227 = x222+(ac/2+2/3*tL22)*cos(pi/180*Inc3);y227 = y222+(ac/2+2/3*tL22)*sin(pi/180*Inc3); x232 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+2*h22+h23)*sin(Inc3*pi/180); y232 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+2*h22+h23)*cos(Inc3*pi/180); x231 = x232-ac/4*cos(pi/180*Inc3);y231 = y232-ac/4*sin(pi/180*Inc3); x233 = x232+ac/4*cos(pi/180*Inc3);y233 = y232+ac/4*sin(pi/180*Inc3); x234 = x232-(ac/2+2/3*tL23)*cos(pi/180*Inc3);y234 = y232-(ac/2+2/3*tL23)*sin(pi/180*Inc3); x235 = x232-(ac/2+1/3*tL23)*cos(pi/180*Inc3);y235 = y232-(ac/2+1/3*tL23)*sin(pi/180*Inc3); x236 = x232+(ac/2+1/3*tL23)*cos(pi/180*Inc3);y236 = y232+(ac/2+1/3*tL23)*sin(pi/180*Inc3); x237 = x232+(ac/2+2/3*tL23)*cos(pi/180*Inc3);y237 = y232+(ac/2+2/3*tL23)*sin(pi/180*Inc3); x242 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+2*h22+2*h23+h24)*sin(Inc3*pi/180); y242 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+2*h22+2*h23+h24)*cos(Inc3*pi/180); x241 = x242-ac/4*cos(pi/180*Inc3);y241 = y242-ac/4*sin(pi/180*Inc3); x243 = x242+ac/4*cos(pi/180*Inc3);y243 = y242+ac/4*sin(pi/180*Inc3); x244 = x242-(ac/2+2/3*tL24)*cos(pi/180*Inc3);y244 = y242-(ac/2+2/3*tL24)*sin(pi/180*Inc3); x245 = x242-(ac/2+1/3*tL24)*cos(pi/180*Inc3);y245 = y242-(ac/2+1/3*tL24)*sin(pi/180*Inc3); x246 = x242+(ac/2+1/3*tL24)*cos(pi/180*Inc3);y246 = y242+(ac/2+1/3*tL24)*sin(pi/180*Inc3); x247 = x242+(ac/2+2/3*tL24)*cos(pi/180*Inc3);y247 = y242+(ac/2+2/3*tL24)*sin(pi/180*Inc3); x252 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+h25)*sin(Inc3*pi/180); y252 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+h25)*cos(Inc3*pi/180); x251 = x252-ac/4*cos(pi/180*Inc3);y251 = y252-ac/4*sin(pi/180*Inc3); x253 = x252+ac/4*cos(pi/180*Inc3);y253 = y252+ac/4*sin(pi/180*Inc3); x254 = x252-(ac/2+2/3*tL25)*cos(pi/180*Inc3);y254 = y252-(ac/2+2/3*tL25)*sin(pi/180*Inc3); x255 = x252-(ac/2+1/3*tL25)*cos(pi/180*Inc3);y255 = y252-(ac/2+1/3*tL25)*sin(pi/180*Inc3); x256 = x252+(ac/2+1/3*tL25)*cos(pi/180*Inc3);y256 = y252+(ac/2+1/3*tL25)*sin(pi/180*Inc3); x257 = x252+(ac/2+2/3*tL25)*cos(pi/180*Inc3);y257 = y252+(ac/2+2/3*tL25)*sin(pi/180*Inc3); x262 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+2*h25+h26)*sin(Inc3*pi/180); y262 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+2*h25+h26)*cos(Inc3*pi/180); x261 = x262-ac/4*cos(pi/180*Inc3);y261 = y262-ac/4*sin(pi/180*Inc3); x263 = x262+ac/4*cos(pi/180*Inc3);y263 = y262+ac/4*sin(pi/180*Inc3); x264 = x262-(ac/2+2/3*tL26)*cos(pi/180*Inc3);y264 = y262-(ac/2+2/3*tL26)*sin(pi/180*Inc3); x265 = x262-(ac/2+1/3*tL26)*cos(pi/180*Inc3);y265 = y262-(ac/2+1/3*tL26)*sin(pi/180*Inc3); x266 = x262+(ac/2+1/3*tL26)*cos(pi/180*Inc3);y266 = y262+(ac/2+1/3*tL26)*sin(pi/180*Inc3); x267 = x262+(ac/2+2/3*tL26)*cos(pi/180*Inc3);y267 = y262+(ac/2+2/3*tL26)*sin(pi/180*Inc3); x272 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+2*h25+2*h26+h27)*sin(Inc3*pi/180); y272 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+2*h25+2*h26+h27)*cos(Inc3*pi/180); x271 = x272-ac/4*cos(pi/180*Inc3);y271 = y272-ac/4*sin(pi/180*Inc3); x273 = x272+ac/4*cos(pi/180*Inc3);y273 = y272+ac/4*sin(pi/180*Inc3); x274 = x272-(ac/2+2/3*tL27)*cos(pi/180*Inc3);y274 = y272-(ac/2+2/3*tL27)*sin(pi/180*Inc3); x275 = x272-(ac/2+1/3*tL27)*cos(pi/180*Inc3);y275 = y272-(ac/2+1/3*tL27)*sin(pi/180*Inc3); x276 = x272+(ac/2+1/3*tL27)*cos(pi/180*Inc3);y276 = y272+(ac/2+1/3*tL27)*sin(pi/180*Inc3); x277 = x272+(ac/2+2/3*tL27)*cos(pi/180*Inc3);y277 = y272+(ac/2+2/3*tL27)*sin(pi/180*Inc3); x282 = kLat*cos(alpL)+Rcc-(Rcc+be/2)*cos(Inc3*pi/180)+0.5*(h19+ce)*sin(Inc3*pi/180)+0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+2*h25+2*h26+2*h27+h28)*sin(Inc3*pi/180); y282 = -(ca/2+cb+cc+cf)-(Rcc+be/2)*sin(Inc3*pi/180)-kLat*sin(alpL)-0.5*(h19+ce)*cos(Inc3*pi/180)-0.5*(h19+2*h20+2*h21+2*h22+2*h23+2*h24+2*h25+2*h26+2*h27+h28)*cos(Inc3*pi/180); x281 = x282-ac/4*cos(pi/180*Inc3);y281 = y282-ac/4*sin(pi/180*Inc3); x283 = x282+ac/4*cos(pi/180*Inc3);y283 = y282+ac/4*sin(pi/180*Inc3); x284 = x282-(ac/2+2/3*tL28)*cos(pi/180*Inc3);y284 = y282-(ac/2+2/3*tL28)*sin(pi/180*Inc3); x285 = x282-(ac/2+1/3*tL28)*cos(pi/180*Inc3);y285 = y282-(ac/2+1/3*tL28)*sin(pi/180*Inc3); x286 = x282+(ac/2+1/3*tL28)*cos(pi/180*Inc3);y286 = y282+(ac/2+1/3*tL28)*sin(pi/180*Inc3); x287 = x282+(ac/2+2/3*tL28)*cos(pi/180*Inc3);y287 = y282+(ac/2+2/3*tL28)*sin(pi/180*Inc3); %Midline Coordinates Matrices coord1a = [x11;y1];coord1b=[x12;y1];coord1c=[x13;y1];coord1d=[x14;y1];coord1e=[x15;y1];coord1f=[x16;y1];coord1g=[x17;y1]; coord2a = [x21;y2];coord2b=[x22;y2];coord2c=[x23;y2];coord2d=[x24;y2];coord2e=[x25;y2];coord2f=[x26;y2];coord2g=[x27;y2]; coord3a = [x31;y3];coord3b=[x32;y3];coord3c=[x33;y3];coord3d=[x34;y3];coord3e=[x35;y3];coord3f=[x36;y3];coord3g=[x37;y3]; coord4a = [x41;y4];coord4b=[x42;y4];coord4c=[x43;y4];coord4d=[x44;y4];coord4e=[x45;y4];coord4f=[x46;y4];coord4g=[x47;y4]; coord5a = [x51;y5];coord5b=[x52;y5];coord5c=[x53;y5];coord5d=[x54;y5];coord5e=[x55;y5];coord5f=[x56;y5];coord5g=[x57;y5]; coord6a = [x61;y6];coord6b=[x62;y6];coord6c=[x63;y6];coord6d=[x64;y6];coord6e=[x65;y6];coord6f=[x66;y6];coord6g=[x67;y6]; coord7a = [x71;y7];coord7b=[x72;y7];coord7c=[x73;y7];coord7d=[x74;y7];coord7e=[x75;y7];coord7f=[x76;y7];coord7g=[x77;y7]; coord8a = [x81;y8];coord8b=[x82;y8];coord8c=[x83;y8];coord8d=[x84;y8];coord8e=[x85;y8];coord8f=[x86;y8];coord8g=[x87;y8]; coord9a = [x91;y9];coord9b=[x92;y9];coord9c=[x93;y9];coord9d=[x94;y9];coord9e=[x95;y9];coord9f=[x96;y9];coord9g=[x97;y9]; coord10a = [x101;y10];coord10b=[x102;y10];coord10c=[x103;y10];coord10d=[x104;y10];coord10e=[x105;y10];coord10f=[x106;y10];coord10g=[x107;y10]; coord11a = [x111;y11];coord11b=[x112;y11];coord11c=[x113;y11];coord11d=[x114;y11];coord11e=[x115;y11];coord11f=[x116;y11];coord11g=[x117;y11]; coord12a = [x121;y12];coord12b=[x122;y12];coord12c=[x123;y12];coord12d=[x124;y12];coord12e=[x125;y12];coord12f=[x126;y12];coord12g=[x127;y12]; coord13a = [x131;y13];coord13b=[x132;y13];coord13c=[x133;y13];coord13d=[x134;y13];coord13e=[x135;y13];coord13f=[x136;y13];coord13g=[x137;y13]; coord14a = [x141;y14];coord14b=[x142;y14];coord14c=[x143;y14];coord14d=[x144;y14];coord14e=[x145;y14];coord14f=[x146;y14];coord14g=[x147;y14]; coord15a = [x151;y15];coord15b=[x152;y15];coord15c=[x153;y15];coord15d=[x154;y15];coord15e=[x155;y15];coord15f=[x156;y15];coord15g=[x157;y15]; coord16a = [x161;y161];coord16b=[x162;y162];coord16c=[x163;y163];coord16d=[x164;y164];coord16e=[x165;y165];coord16f=[x166;y166];coord16g=[x167;y167]; coord17a = [x171;y171];coord17b=[x172;y172];coord17c=[x173;y173];coord17d=[x174;y174];coord17e=[x175;y175];coord17f=[x176;y176];coord17g=[x177;y177]; coord18a = [x181;y181];coord18b=[x182;y182];coord18c=[x183;y183];coord18d=[x184;y184];coord18e=[x185;y185];coord18f=[x186;y186];coord18g=[x187;y187]; coord19a = [x191;y191];coord19b=[x192;y192];coord19c=[x193;y193];coord19d=[x194;y194];coord19e=[x195;y195];coord19f=[x196;y196];coord19g=[x197;y197]; coord20a = [x201;y201];coord20b=[x202;y202];coord20c=[x203;y203];coord20d=[x204;y204];coord20e=[x205;y205];coord20f=[x206;y206];coord20g=[x207;y207]; coord21a = [x211;y211];coord21b=[x212;y212];coord21c=[x213;y213];coord21d=[x214;y214];coord21e=[x215;y215];coord21f=[x216;y216];coord21g=[x217;y217]; coord22a = [x221;y221];coord22b=[x222;y222];coord22c=[x223;y223];coord22d=[x224;y224];coord22e=[x225;y225];coord22f=[x226;y226];coord22g=[x227;y227]; coord23a = [x231;y231];coord23b=[x232;y232];coord23c=[x233;y233];coord23d=[x234;y234];coord23e=[x235;y235];coord23f=[x236;y236];coord23g=[x237;y237]; coord24a = [x241;y241];coord24b=[x242;y242];coord24c=[x243;y243];coord24d=[x244;y244];coord24e=[x245;y245];coord24f=[x246;y246];coord24g=[x247;y247]; coord25a = [x251;y251];coord25b=[x252;y252];coord25c=[x253;y253];coord25d=[x254;y254];coord25e=[x255;y255];coord25f=[x256;y256];coord25g=[x257;y257]; coord26a = [x261;y261];coord26b=[x262;y262];coord26c=[x263;y263];coord26d=[x264;y264];coord26e=[x265;y265];coord26f=[x266;y266];coord26g=[x267;y267]; coord27a = [x271;y271];coord27b=[x272;y272];coord27c=[x273;y273];coord27d=[x274;y274];coord27e=[x275;y275];coord27f=[x276;y276];coord27g=[x277;y277]; coord28a = [x281;y281];coord28b=[x282;y282];coord28c=[x283;y283];coord28d=[x284;y284];coord28e=[x285;y285];coord28f=[x286;y286];coord28g=[x287;y287]; Tt=30; %total running time in days %average velocity at the midpoint of each layer..livelink with comsol in %m/s for k=1:Tt vfmid1(k) = 1/7*(mphinterp(model,'u','coord',coord1a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord1b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord1c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord1d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord1e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord1f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord1g,'dataset','dset2','t',k)); vfmid2(k) = 1/7*(mphinterp(model,'u','coord',coord2a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord2b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord2c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord2d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord2e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord2f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord2g,'dataset','dset2','t',k)); vfmid3(k) = 1/7*(mphinterp(model,'u','coord',coord3a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord3b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord3c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord3d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord3e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord3f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord3g,'dataset','dset2','t',k)); vfmid4(k) = 1/7*(mphinterp(model,'u','coord',coord4a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord4b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord4c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord4d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord4e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord4f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord4g,'dataset','dset2','t',k)); vfmid5(k) = 1/7*(mphinterp(model,'u','coord',coord5a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord5b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord5c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord5d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord5e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord5f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord5g,'dataset','dset2','t',k)); vfmid6(k) = 1/7*(mphinterp(model,'u','coord',coord6a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord6b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord6c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord6d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord6e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord6f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord6g,'dataset','dset2','t',k)); vfmid7(k) = 1/7*(mphinterp(model,'u','coord',coord7a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord7b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord7c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord7d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord7e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord7f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord7g,'dataset','dset2','t',k)); vfmid8(k) = 1/7*(mphinterp(model,'u','coord',coord8a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord8b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord8c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord8d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord8e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord8f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord8g,'dataset','dset2','t',k)); vfmid9(k) = 1/7*(mphinterp(model,'u','coord',coord9a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord9b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord9c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord9d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord9e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord9f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord9g,'dataset','dset2','t',k)); vfmid10(k) = 1/7*(mphinterp(model,'u','coord',coord10a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord10b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord10c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord10d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord10e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord10f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord10g,'dataset','dset2','t',k)); vfmid11(k) = 1/7*(mphinterp(model,'u','coord',coord11a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord11b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord11c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord11d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord11e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord11f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord11g,'dataset','dset2','t',k)); vfmid12(k) = 1/7*(mphinterp(model,'u','coord',coord12a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord12b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord12c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord12d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord12e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord12f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord12g,'dataset','dset2','t',k)); vfmid13(k) = 1/7*(mphinterp(model,'u','coord',coord13a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord13b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord13c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord13d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord13e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord13f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord13g,'dataset','dset2','t',k)); vfmid14(k) = 1/7*(mphinterp(model,'u','coord',coord14a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord14b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord14c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord14d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord14e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord14f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord14g,'dataset','dset2','t',k)); vfmid15(k) = 1/7*(mphinterp(model,'u','coord',coord15a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord15b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord15c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord15d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord15e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord15f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord15g,'dataset','dset2','t',k)); vfmid16(k) = 1/7*(mphinterp(model,'u','coord',coord16a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord16b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord16c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord16d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord16e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord16f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord16g,'dataset','dset2','t',k)); vfmid17(k) = 1/7*(mphinterp(model,'u','coord',coord17a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord17b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord17c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord17d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord17e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord17f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord17g,'dataset','dset2','t',k)); vfmid18(k) = 1/7*(mphinterp(model,'u','coord',coord18a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord18b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord18c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord18d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord18e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord18f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord18g,'dataset','dset2','t',k)); vfmid19(k) = 1/7*(mphinterp(model,'u','coord',coord19a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord19b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord19c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord19d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord19e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord19f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord19g,'dataset','dset2','t',k)); vfmid20(k) = 1/7*(mphinterp(model,'u','coord',coord20a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord20b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord20c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord20d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord20e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord20f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord20g,'dataset','dset2','t',k)); vfmid21(k) = 1/7*(mphinterp(model,'u','coord',coord21a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord21b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord21c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord21d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord21e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord21f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord21g,'dataset','dset2','t',k)); vfmid22(k) = 1/7*(mphinterp(model,'u','coord',coord22a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord22b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord22c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord22d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord22e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord22f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord22g,'dataset','dset2','t',k)); vfmid23(k) = 1/7*(mphinterp(model,'u','coord',coord23a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord23b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord23c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord23d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord23e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord23f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord23g,'dataset','dset2','t',k)); vfmid24(k) = 1/7*(mphinterp(model,'u','coord',coord24a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord24b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord24c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord24d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord24e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord24f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord24g,'dataset','dset2','t',k)); vfmid25(k) = 1/7*(mphinterp(model,'u','coord',coord25a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord25b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord25c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord25d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord25e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord25f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord25g,'dataset','dset2','t',k)); vfmid26(k) = 1/7*(mphinterp(model,'u','coord',coord26a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord26b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord26c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord26d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord26e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord26f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord26g,'dataset','dset2','t',k)); vfmid27(k) = 1/7*(mphinterp(model,'u','coord',coord27a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord27b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord27c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord27d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord27e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord27f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord27g,'dataset','dset2','t',k)); vfmid28(k) = 1/7*(mphinterp(model,'u','coord',coord28a,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord28b,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord28c,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord28d,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord28e,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord28f,'dataset','dset2','t',k)+mphinterp(model,'u','coord',coord28g,'dataset','dset2','t',k)); end %average velocities at midpoints in matrix in m/s vfmid=[vfmid1' vfmid2' vfmid3' vfmid4' vfmid5' vfmid6' vfmid7' vfmid8' vfmid9' vfmid10' vfmid11' vfmid12' vfmid13' vfmid14'... vfmid15' vfmid16' vfmid17' vfmid18' vfmid19' vfmid20' vfmid21' vfmid22' vfmid23' vfmid24' vfmid25' vfmid26' vfmid27' vfmid28']; %average pressure at the midpoint of each layer...livelink with comsol in %N/m^2 for k=1:Tt pwmid1(k) = 1/7*(mphinterp(model,'p','coord',coord1a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord1b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord1c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord1d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord1e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord1f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord1g,'dataset','dset2','t',k)); pwmid2(k) = 1/7*(mphinterp(model,'p','coord',coord2a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord2b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord2c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord2d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord2e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord2f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord2g,'dataset','dset2','t',k)); pwmid3(k) = 1/7*(mphinterp(model,'p','coord',coord3a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord3b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord3c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord3d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord3e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord3f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord3g,'dataset','dset2','t',k)); pwmid4(k) = 1/7*(mphinterp(model,'p','coord',coord4a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord4b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord4c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord4d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord4e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord4f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord4g,'dataset','dset2','t',k)); pwmid5(k) = 1/7*(mphinterp(model,'p','coord',coord5a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord5b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord5c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord5d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord5e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord5f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord5g,'dataset','dset2','t',k)); pwmid6(k) = 1/7*(mphinterp(model,'p','coord',coord6a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord6b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord6c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord6d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord6e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord6f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord6g,'dataset','dset2','t',k)); pwmid7(k) = 1/7*(mphinterp(model,'p','coord',coord7a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord7b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord7c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord7d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord7e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord7f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord7g,'dataset','dset2','t',k)); pwmid8(k) = 1/7*(mphinterp(model,'p','coord',coord8a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord8b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord8c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord8d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord8e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord8f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord8g,'dataset','dset2','t',k)); pwmid9(k) = 1/7*(mphinterp(model,'p','coord',coord9a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord9b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord9c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord9d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord9e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord9f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord9g,'dataset','dset2','t',k)); pwmid10(k) = 1/7*(mphinterp(model,'p','coord',coord10a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord10b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord10c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord10d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord10e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord10f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord10g,'dataset','dset2','t',k)); pwmid11(k) = 1/7*(mphinterp(model,'p','coord',coord11a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord11b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord11c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord11d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord11e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord11f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord11g,'dataset','dset2','t',k)); pwmid12(k) = 1/7*(mphinterp(model,'p','coord',coord12a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord12b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord12c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord12d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord12e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord12f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord12g,'dataset','dset2','t',k)); pwmid13(k) = 1/7*(mphinterp(model,'p','coord',coord13a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord13b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord13c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord13d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord13e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord13f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord13g,'dataset','dset2','t',k)); pwmid14(k) = 1/7*(mphinterp(model,'p','coord',coord14a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord14b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord14c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord14d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord14e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord14f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord14g,'dataset','dset2','t',k)); pwmid15(k) = 1/7*(mphinterp(model,'p','coord',coord15a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord15b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord15c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord15d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord15e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord15f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord15g,'dataset','dset2','t',k)); pwmid16(k) = 1/7*(mphinterp(model,'p','coord',coord16a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord16b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord16c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord16d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord16e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord16f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord16g,'dataset','dset2','t',k)); pwmid17(k) = 1/7*(mphinterp(model,'p','coord',coord17a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord17b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord17c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord17d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord17e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord17f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord17g,'dataset','dset2','t',k)); pwmid18(k) = 1/7*(mphinterp(model,'p','coord',coord18a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord18b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord18c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord18d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord18e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord18f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord18g,'dataset','dset2','t',k)); pwmid19(k) = 1/7*(mphinterp(model,'p','coord',coord19a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord19b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord19c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord19d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord19e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord19f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord19g,'dataset','dset2','t',k)); pwmid20(k) = 1/7*(mphinterp(model,'p','coord',coord20a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord20b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord20c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord20d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord20e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord20f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord20g,'dataset','dset2','t',k)); pwmid21(k) = 1/7*(mphinterp(model,'p','coord',coord21a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord21b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord21c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord21d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord21e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord21f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord21g,'dataset','dset2','t',k)); pwmid22(k) = 1/7*(mphinterp(model,'p','coord',coord22a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord22b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord22c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord22d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord22e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord22f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord22g,'dataset','dset2','t',k)); pwmid23(k) = 1/7*(mphinterp(model,'p','coord',coord23a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord23b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord23c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord23d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord23e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord23f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord23g,'dataset','dset2','t',k)); pwmid24(k) = 1/7*(mphinterp(model,'p','coord',coord24a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord24b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord24c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord24d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord24e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord24f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord24g,'dataset','dset2','t',k)); pwmid25(k) = 1/7*(mphinterp(model,'p','coord',coord25a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord25b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord25c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord25d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord25e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord25f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord25g,'dataset','dset2','t',k)); pwmid26(k) = 1/7*(mphinterp(model,'p','coord',coord26a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord26b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord26c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord26d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord26e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord26f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord26g,'dataset','dset2','t',k)); pwmid27(k) = 1/7*(mphinterp(model,'p','coord',coord27a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord27b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord27c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord27d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord27e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord27f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord27g,'dataset','dset2','t',k)); pwmid28(k) = 1/7*(mphinterp(model,'p','coord',coord28a,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord28b,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord28c,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord28d,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord28e,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord28f,'dataset','dset2','t',k)+mphinterp(model,'p','coord',coord28g,'dataset','dset2','t',k)); end %average pressures at midpoints in matrix ...in Psi pwmid=0.000145038*[pwmid1' pwmid2' pwmid3' pwmid4' pwmid5' pwmid6' pwmid7' pwmid8' pwmid9' pwmid10' pwmid11' pwmid12' pwmid13' pwmid14'... pwmid15' pwmid16' pwmid17' pwmid18' pwmid19' pwmid20' pwmid21' pwmid22' pwmid23' pwmid24' pwmid25' pwmid26' pwmid27' pwmid28']; pwmidp=pwmid; %wellbore flowing pressure at producing zones % MATLAB INPUT DATA nL = 15; % Total Number of Sections (openhole and non-openhole) nnL = 28; % Number of Openhole Sections nc=4; %non-openhole sections open to wellbore flow nst=1; %Beginning of Openhole layer nend=15;%End of Openhole layer C0 = 100*[9.71 3.041 4.359 3.546 2.580 1.699 3.138 2.644 1.291 1.657 9.71 3.041 4.359 3.546 2.580 1.699 3.138 2.644 1.291 1.657... 9.71 3.041 4.359 9.71 3.041 4.359 3.546 2.580 1.699 3.138 2.644 1.291 1.657]; % Cohesive Strength of the rock Phi = 0.5*[0.6624 0.6306 0.569313 0.5693 0.611 0.611 0.611 0.611 0.611 0.611 0.6624 0.6306 0.569313 0.5693 0.611 0.611 0.611 0.611 0.611 0.611... 0.6624 0.6306 0.569313 0.6624 0.6306 0.569313 0.5693 0.611 0.611 0.611 0.611 0.611 0.611]; % Internal friction angle of the rocks mu = tan(Phi); % coefficient of internal friction P0 = 0.5*[2820 2828 2834 2840 2849 2851 2863 2864 2866 2878 2820 2828 2834 2840 2849 2851 2863 2864 2866 2878 ... 2820 2828 2834 2840 2849 2851 2863 2864 2866 2878 2864 2866 2878]; % Pore pressure of the formations v = 0.3*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % Poisson Ratio %a = 6*ones(1,nnL); %Initial Wellbore Radius before breakout por = 0.1*[5 5 5 5 4 4 4 3 3 3 3 2 2 2 2 1.5 1.5 1.5 1.4 1.3 1.3 1.25 1.2 ... 1.25 1.3 1.13 1.2 1.1 1.25 1.3 1.13 1.2 1.1]; %porosity of each openhole layers (dimensionless) porc=0.1*[2 1.5 3 0.1];%Blowhole porosity in the cement; Noting that the Conductor Casing's is the first element and Liner is the last (in that order) KICf=1000*ones(1,nL);%Undamaged fracture toughness of the formations; from top to bottom KICc=3000*ones(1,nc);%Undamaged fracture toughness of cement...from liner-cement to conductor casing-cement %Damaged Fracture toughness of cement...noting that nst and nend are not %applicable here % hpf =12*[hc(1) hc(2) hc(3) hc(4) h1 h2 h3 h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16 h17 h18 h19 h20... % h21 h22 h23 h24 h25 h26 h27 h28]; %in inches % h = 40*ones(1,nL); for j=1:nc KICcd(j)=KICc(j).*(1-porc(j)).^(3/2); end KICs=1e10*ones(1,nc);%Fracture Toughness of casings...assume a very high number %Insitu Stresses in psi SV = 2000*[5.778 5.795 5.808 5.808 5.822 5.843 5.848 5.875 5.876 5.883 5.778 5.795 5.808 5.808 5.822 5.843 5.848 5.875 5.876 5.883... 5.778 5.795 5.808 5.778 5.795 5.808 5.808 5.822 5.843 5.848 5.875 5.876 5.883]; %Overburden SH = 2000*[6.588 6.608 6.623 6.624 6.640 6.664 6.67 6.702 6.703 6.710 6.588 6.608 6.623 6.624 6.640 6.664 6.67 6.702 6.703 6.710... 6.588 6.608 6.623 6.588 6.608 6.623 6.624 6.640 6.664 6.67 6.702 6.703 6.710]; Shm = 2000*[4.088 4.099 4.108 4.108 4.117 4.13 4.134 4.152 4.152 4.157 4.088 4.099 4.108 4.108 4.117 4.13 4.134 4.152 4.152 4.157... 4.088 4.099 4.108 4.088 4.099 4.108 4.108 4.117 4.13 4.134 4.152 4.152 4.157]; %Minimum Horizontal Stress in psi %Well Angles incl = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inchlf1 Inchlf2 Inchlf3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3 pi/180*Inc3]; %inclination angles azi = pi/180*[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 135 135 135 135 135 135 135 135 135 135]; %azimumthal angles % Reservoir(s) Properties kh = [3 3 3 3 2 2 2 1 1 1 1 0.5 0.5 0.5 0.5 0.25 0.25 0.25 0.2 0.2 0.2 0.15 0.15 ... 0.2 0.25 0.11 0.15 0.1]; %in-plane permeability of the openhole sections in mD kv = [3 3 3 3 2 2 2 1 1 1 1 0.5 0.5 0.5 0.5 0.25 0.25 0.25 0.2 0.2 0.2 0.15 0.15 ... 0.2 0.25 0.11 0.15 0.1]; %out of plane permeability of the producing formations in mD pint = P0; %Initial reservoirs pressures in psia mur = [3 3 3 3 2 2 2 1 1 1 1 0.5 0.5 0.5 0.5 0.25 0.25 0.25 0.2 0.2 0.2 0.15 0.15 ... 0.2 0.25 0.11 0.15 0.1]; % Viscosity of reservoirs fluids in centipoise index=1; %Choose 1 if producing formation is sandstones %Choose any other number if producing formation is limestone M = [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ]; %Formation fluid type: 1-gas, 2- liquid inclp=incl;% wellbore inclination angle along the producing formations G = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 2]; %Boundary condition; applicable to horizontal well % 1=noflow boundary conditions up and down, 2=constant pressure % up and noflow boundary down gg = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %Gas specific gravity T = [600 600 600 600 600 600 600 600 600 650 650 650 650 650 650 650 700 700 700 700 700 700 700 750 750 ... 750 750 750];%Formation Temperature in Rankine psc = 14.7; %standard pressure in psia tsc = 520; %standard temperature in Rankine nnp=nend-nst+1; %number of producing openhole intervals z0=ones(1,nnp); %initial z-factors accz = 0.00001*ones(1,nnp);%accuracy for estimating z-factor %Reservoir Dimensions in ft aw = [100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 ... 100 100 100 100 100 100 100 100]; %Reservoir width(can not be zero) bL = [100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 ... 100 100 100 100 100 100 100 100]; %Reservoir Length(can not be zero) %Well Positions in ft x1 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %Initial Longitudinal positions x2 = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %Final Longitudinal Positions zw = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%Vertical distance from the center of horizontal well to the bottom of the formation yw = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];%Lateral distance of the well position to the origin (heel) L = h; %Well Length in each formation %thickness of well penetrated section for j=nst:nend hw(j) = L(j).*cos(incl(j)); end %Estimating Fluid Properties: Viscosity and Density % Stokes Settling Velocity for a single particle; assuming all % particles have same diameter or use mean diameter of the particles dp=1; %mean particle diameter in cm mufw= (2)*0.001; %reservoir fluid initial viscosity (viscosity calculator in centipoise)in Pa.s rhof= 1; % density of reservoir fluid in g/cm^3; rhos=2.3; %mean density of solid ingress to the wellbore in g/cm^3 %CALCULATIONS %Insitu Stresses along Wellbore Orientations for j=nst:nend SX(j)=((SH(j)-P0(j)).*(cos(azi(j))).^2+(Shm(j)-P0(j)).*(sin(azi(j))).^2).*(cos(incl(j))).^2+(SV(j)-P0(j)).*(sin(incl(j))).^2; SY(j) = ((SH(j)-P0(j)).*(sin(azi(j))).^2+(Shm(j)-P0(j)).*(cos(azi(j))).^2); SZZ(j)=((SH(j)-P0(j)).*(cos(azi(j))).^2+(Shm(j)-P0(j)).*(sin(azi(j))).^2).*(sin(incl(j))).^2+(SV(j)-P0(j)).*(cos(incl(j))).^2; SXY(j)=0.5*(Shm(j)-SH(j)).*sin(2*azi(j)).*cos(incl(j)); SXZ(j)=0.5*((SH(j)-P0(j)).*(cos(azi(j))).^2+(Shm(j)-P0(j)).*(sin(azi(j))).^2-(SV(j)-P0(j))).*sin(2*incl(j)); SYZ(j)=0.5*(Shm(j)-SH(j)).*sin(2*azi(j)).*sin(incl(j)); end for j=nst:nend if SX(j)>SY(j) SXv(j)=SX(j); SYv(j)=SY(j); else SXv(j)=SY(j); SYv(j)=SX(j); end end % WELLBORE COLLAPSE PRESSURE FUNCTION FOR EACH LAYER % A.Circular Wellbore % Using Mohr Coulomb Failure Criterion syms thbk pw for k=1:Tt for j=nst:1:nend % Disregarding the out of plane shear stresses; Breakout Initiation Angle thbkfun(j)=0.5*atan(2*SXY(j)./(SXv(j)-SYv(j))); if SX(j)>SY(j) && SXY(j)<=0 thbkso(k,j)=pi/2-thbkfun(j); elseif SX(j)>SY(j) && SXY(j)>=0 thbkso(k,j)=pi/2+thbkfun(j); elseif SY(j)>SX(j) && SXY(j)<=0 thbkso(k,j)=pi-thbkfun(j); elseif SY(j)>SX(j) && SXY(j)>=0 thbkso(k,j)=thbkfun(j); end %Stresses around the wellbore wall at maximum compressive direction SRwbk(k,j)=pwmid(k,j)-P0(j); STwbk(k,j)=(SXv(j)+SYv(j))-2*(SXv(j)-SYv(j)).*cos(2*thbkso(k,j))-4*SXY(j)*sin(2*thbkso(k,j))-pwmid(k,j)+P0(j); SZwbk(k,j)=SZZ(j)-v(j).*(2*(SXv(j)-SYv(j))*cos(2*thbkso(k,j))+4*SXY(j)*sin(2*thbkso(k,j))); %Stresses around the wellbore wall at minimum compressive direction SRwbkm(k,j)=pwmid(k,j)-P0(j); STwbkm(k,j)=(SXv(j)+SYv(j))-2*(SXv(j)-SYv(j)).*cos(2*thbkso(k,j)+pi)-4*SXY(j)*sin(2*thbkso(k,j)+pi)-pwmid(k,j)+P0(j); SZwbkm(k,j)=SZZ(j)-v(j).*(2*(SXv(j)-SYv(j))*cos(2*thbkso(k,j)+pi)+4*SXY(j)*sin(2*thbkso(k,j)+pi)); % Collapse Pressure if (STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j)) %Type 1 Breakout SRw(j)=pw-P0(j); STw(j)=(SXv(j)+SYv(j))-2*(SXv(j)-SYv(j)).*cos(2*thbk)-4*SXY(j)*sin(2*thbk)-pw+P0(j); SZw(j)=SZZ(j)-v(j).*(2*(SXv(j)-SYv(j))*cos(2*thbk)+4*SXY(j)*sin(2*thbk)); STZw(j)=2*(-SXZ(j)*sin(thbk)+SYZ(j)*cos(thbk)); %Considering Out-of-Plane Shear Stresses sigma1w(j)=0.5*(STw(j)+SZw(j)+sqrt((STw(j)-SZw(j)).^2+4*STZw(j).^2)); sigma3w(j)= SRw(j); equation1(j)=0.5*(sigma1w(j)-sigma3w(j)).*cos(Phi(j)); equation2(j)=C0(j)+(0.5*(sigma1w(j)+sigma3w(j))-0.5*(sigma1w(j)-sigma3w(j)).*sin(Phi(j))).*tan(Phi(j)); pwfun{j}=solve(equation1(j)==equation2(j),pw); pwfuns1(j)=pwfun{j}(1); pwfuns2(j)=pwfun{j}(2); %Breakout Direction when Considering Out-of-Plane Shear Stresses thbkso1(k,j)=double(vpasolve(diff(pwfuns1(j),thbk)==0,thbk)); thbkso2(k,j)=double(vpasolve(diff(pwfuns2(j),thbk)==0,thbk)); % Collapse Pressure pwc1(k,j)=double(subs(pwfuns1(j),thbk,thbkso1(k,j))); pwc2(k,j)=double(subs(pwfuns2(j),thbk,thbkso2(k,j))); pwc(k,j)=(0.5*(2*P0(j).*sin(Phi(j))-2*C0(j).*cos(Phi(j))+(2*cos(2*thbkso(k,j))-1).*SXv(j).*(sin(Phi(j))-1)... +4*sin(2*thbkso(k,j)).*SXY(j).*sin(Phi(j))-4*sin(2*thbkso(k,j)).*SXY(j)-(2*cos(2*thbkso(k,j))+1).*SYv(j).*(sin(Phi(j))-1))); elseif (SZwbk(k,j)>STwbk(k,j) && STwbk(k,j)>SRwbk(k,j)) % Type 3 Breakout %Collapse Pressure pwc(k,j)=(-C0(j).*cos(Phi(j))+sin(Phi(j)).*(P0(j)-SZZ(j))+P0(j)+2*v(j).*(sin(Phi(j))-1).*cos(2*thbkso(k,j)).(SXv(j)-SYv(j))+2*SXY(j).*sin(2*thbkso(k,j))+SZZ(j))./(1+sin(Phi(j))); end STZwc(k,j)=2*(-SXZ(j)*sin(thbkso(k,j))+SYZ(j)*cos(thbkso(k,j))); end end % CHECK FOR COLLAPSE IN EACH LAYER AND ESTIMATE BREAKOUT DIMENSIONS %Estimating Breakout Depth in the direction of maximum compression(rb) syms rD rDn for k=1:Tt for j=nst:1:nend if pwmid(k,j)<=pwc(k,j) SRR(k,j)=0.5*(SXv(j)+SYv(j))*(1-rD^2)+0.5*(SXv(j)-SYv(j)).*(1+(3*rD^4)-4*rD^2).*cos(2*thbkso(k,j))... +SXY(j).*(1+(3*rD^4)-4*rD^2).*sin(2*thbkso(k,j))+(rD^2)*(pwmid(k,j)-P0(j)); STT(k,j)=0.5*(SXv(j)+SYv(j)).*(1+rD^2)-0.5*(SXv(j)-SYv(j)).*(1+3*rD^4).*cos(2*thbkso(k,j))... -SXY(j).*(1+3*rD^4).*sin(2*thbkso(k,j))-(rD^2).*(pwmid(k,j)-P0(j)); SRT(k,j)=(0.5*(SXv(j)-SYv(j)).*sin(2*thbkso(k,j))+SXY(j).*cos(2*thbkso(k,j))).*(1-(3*rD^4)+(2*rD^2)); if (STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j)) sigma1A(k,j)=0.5*(STT(k,j)+SRR(k,j))+0.5*((STT(k,j)-SRR(k,j)).^2+4*SRT(k,j).^2).^0.5; sigma3A(k,j)=0.5*(STT(k,j)+SRR(k,j))-0.5*((STT(k,j)-SRR(k,j)).^2+4*SRT(k,j).^2).^0.5; elseif (SZwbk(k,j)>STwbk(k,j) && STwbk(k,j)>SRwbk(k,j)) sigma1A(k,j)=SZZ(j)-v(j).*(2*(SXv(j)-SYv(j))*cos(2*thbkso(k,j))+4*SXY(j)*sin(2*thbkso(k,j)))*rD^2; sigma3A(k,j)=0.5*(STT(k,j)+SRR(k,j))-0.5*((STT(k,j)-SRR(k,j)).^2+4*SRT(k,j).^2).^0.5; end %C0fun(j) = sqrt(1+mu(j).^2).*sqrt(0.25*(STT(j)-SRR(j)).^2+SRT(j).^2)-mu(j).*(0.5*(STT(j)+SRR(j))); C0fun(k,j)=0.5*((sigma1A(k,j).*((1+mu(j).^2).^0.5-mu(j)))-(sigma3A(k,j).*((1+mu(j).^2).^0.5+mu(j)))); %C0fun(j) = 0.5*(-sigma3A(j).*(mu(j)+sqrt(1+mu(j)))+sigma1A(j)./(mu(j)+sqrt(1+mu(j)))); C0funtest=@(rD)eval(C0fun(k,j)-C0(j)); options=optimoptions('fsolve','Display','iter'); rD1(k,j)=real(fsolve(C0funtest,0.7,options)); if rD1(k,j)<0 rbmx1(k,j)=a(j); elseif rD1(k,j)>1 rbmx1(k,j)=a(j); else rbmx1(k,j)= a(j)./rD1(k,j); %breakout depth in the direction of max. compression end else rbmx1(k,j)=a(j); end end end % Estimating the Breakout Width in Max Compresion Direction for k=1:Tt for j=nst:1:nend if pwmid(k,j)<=pwc(j) if (STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j)) sigma3b(k,j)=pwmid(k,j)-P0(j); A(k,j)=(2*C0(j)+sigma3b(k,j).*((1+mu(j).^2).^0.5+mu(j)))./((1+mu(j).^2).^0.5-mu(j)); B(j)=4*SXY(j); C(k,j)=SXv(j)+SYv(j)-A(k,j)-pwmid(k,j)-P0(j); D(j)=2*(SXv(j)-SYv(j)); omega(j)=atan(B(j)./D(j)); lambda(j)=sqrt(D(j).^2+B(j).^2); %Analytical Solution thetab1(k,j)=real(0.5*(acos(C(k,j)./lambda(j))+omega(j))); %idtheta1 = thetab1>=pi/2; %thetab1(idtheta1)=pi/2; %elseif thetab1(j)<0 %thetab1(j)=pi/2; %else %thetab1(j)=pi/2; elseif (SZwbk(k,j)>STwbk(k,j) && STwbk(k,j)>SRwbk(k,j)) C(k,j)=(2*C0(j)+(pwmid(k,j)-P0(j)).*(sqrt(1+mu(j).^2)+mu(j)))./(sqrt(1+mu(j).^2)-mu(j)); D(j)=SZZ(j)-2*v(j).*(SXv(j)-SYv(j)); B(j)=4*v(j).*SXY(j); omega(j)=atan(B(j)./D(j)); lambda(j)=sqrt(D(j).^2+B(j).^2); thetab1(k,j)=real(0.5*(acos(C(k,j)./lambda(j))-omega(j))); end if rbmx1(k,j)>a(j) if thetab1(k,j)<0 thetabkt(k,j)=2*pi+thetab1(k,j); % Breakout width elseif thetab1(k,j)>2*pi thetabkt(k,j)=2*pi; else thetabkt(k,j)=thetab1(k,j); end else thetabkt(k,j)=0; end else thetabkt(k,j)=0; end end end % Estimating the Breakout Depth in the minimum compression direction syms rDc rDmn for k=1:Tt for j=nst:1:nend if pwmid(k,j)<=pwc(j) SRRm(k,j)=0.5*(SXv(j)+SYv(j))*(1-rDc^2)+0.5*(SXv(j)-SYv(j)).*(1+(3*rDc^4)-4*rDc^2).*cos(2*(thbkso(k,j)+1.5708))... +SXY(j).*(1+(3*rDc^4)-4*rDc^2).*sin(2*(thbkso(k,j)+1.5708))+(rDc^2)*(pwmid(k,j)-P0(j)); STTm(k,j)=0.5*(SXv(j)+SYv(j)).*(1+rDc^2)-0.5*(SXv(j)-SYv(j)).*(1+3*rDc^4).*cos(2*(thbkso(k,j)+1.5708))... -SXY(j).*(1+3*rDc^4).*sin(2*(thbkso(k,j)+1.5708))-(rDc^2).*(pwmid(k,j)-P0(j)); SRTm(k,j)=(0.5*(SXv(j)-SYv(j)).*sin(2*(thbkso(k,j)+1.5708))+SXY(j).*cos(2*(thbkso(k,j)+1.5708))).*(1-(3*rDc^4)+(2*rDc^2)); if (STwbkm(k,j)>SZwbkm(k,j) && SZwbkm(k,j)>SRwbkm(k,j)) sigma1Am(k,j)=0.5*(STTm(k,j)+SRRm(k,j))+0.5*((STTm(k,j)-SRRm(k,j)).^2+4*SRTm(k,j).^2).^0.5; sigma3Am(k,j)=0.5*(STTm(k,j)+SRRm(k,j))-0.5*((STTm(k,j)-SRRm(k,j)).^2+4*SRTm(k,j).^2).^0.5; elseif (SZwbkm(k,j)>STwbkm(k,j) && STwbkm(k,j)>SRwbkm(k,j)) sigma1Am(k,j)=SZZ(j)-v(j).*(2*(SXv(j)-SYv(j))*cos(2*thbkso(k,j)+pi)+4*SXY(j)*sin(2*thbkso(k,j)+pi))*rDc^2; sigma3Am(k,j)=0.5*(STTm(k,j)+SRRm(k,j))-0.5*((STTm(k,j)-SRRm(k,j)).^2+4*SRTm(k,j).^2).^0.5; end %C0fun(j) = sqrt(1+mu(j).^2).*sqrt(0.25*(STT(j)-SRR(j)).^2+SRT(j).^2)-mu(j).*(0.5*(STT(j)+SRR(j))); C0funm(k,j)=0.5*((sigma1Am(k,j).*((1+mu(j).^2).^0.5-mu(j)))-(sigma3Am(k,j).*((1+mu(j).^2).^0.5+mu(j)))); %C0funm(j) = 0.5*(-sigma3A(j).*(mu(j)+sqrt(1+mu(j)))+sigma1A(j)./(mu(j)+sqrt(1+mu(j)))); C0funmtest=@(rDc)eval(C0funm(k,j)-C0(j)); options=optimoptions('fsolve','Display','iter'); rD1m(k,j)=real(fsolve(C0funmtest,0.7,options)); if rD1m(k,j)<0 rbmxc(k,j)=a(j); rD1m(k,j)=1; elseif rD1m(k,j)>1 rbmxc(k,j)=a(j); rD1m(k,j)=1; else rbmxc(k,j)= a(j)./rD1m(k,j); %breakout depth in the direction of max. compression end else rbmxc(k,j)=a(j); end end end % Estimating the Breakout Width in Minimum Compresion Direction when Types 1 and 3 Occur Together for k=1:Tt for j=nst:1:nend if pwmid(k,j)<=pwc(j) if thetabkt(k,j)STwbkm(k,j) && STwbkm(k,j)>SRwbkm(k,j)) % pwcm(j)=(-C0(j).*cos(Phi(j))+sin(Phi(j)).*(P0(j)-SZZ(j))+P0(j)+2*v(j).*(sin(Phi(j))-1).*cos(2*thbkso(j)+pi).(SXv(j)-SYv(j))+2*SXY(j).*sin(2*thbkso(j)+pi)+SZZ(j))./(1+sin(Phi(j))); Cm(k,j)=(2*C0(j)+(pwmid(k,j)-P0(j)).*(sqrt(1+mu(j).^2)+mu(j)))./(sqrt(1+mu(j).^2)-mu(j)); Dm(j)=SZZ(j)-2*v(j).*(SXv(j)-SYv(j)); Bm(j)=4*v(j).*SXY(j); omegam(j)=atan(Bm(j)./Dm(j)); lambdam(j)=sqrt(Dm(j).^2+Bm(j).^2); thetab1m(k,j)=real(0.5*(acos(Cm(k,j)./lambdam(j))-omegam(j))); elseif (STwbkm(k,j)>SZwbkm(k,j) && SZwbkm(k,j)>SRwbkm(k,j)) sigma3bm(k,j)=pwmid(k,j)-P0(j); Am(k,j)=(2*C0(j)+sigma3bm(k,j).*((1+mu(j).^2).^0.5+mu(j)))./((1+mu(j).^2).^0.5-mu(j)); Bm(j)=4*SXY(j); Cm(k,j)=SXv(j)+SYv(j)-Am(k,j)-pwmid(k,j)-P0(j); Dm(j)=2*(SXv(j)-SYv(j)); omegam(j)=atan(Bm(j)./Dm(j)); lambdam(j)=sqrt(Dm(j).^2+Bm(j).^2); %Analytical Solution thetab1m(k,j)=real(0.5*(acos(Cm(k,j)./lambdam(j))+omegam(j))); else thetab1m(k,j)=0; end if rbmxc(k,j)>a(j) if thetab1m(k,j)<0 thetabktm(k,j)=2*pi+thetab1m(k,j); % Breakout width elseif thetab1m(k,j)>2*pi thetabktm(k,j)=2*pi; else thetabktm(k,j)=thetab1m(k,j); end else thetabktm(k,j)=0; end end else thetabktm(k,j)=0; end end end % B. Elliptic Wellbore syms thbke thbkn for k=1:Tt for j=nst:1:nend if rbmx1(k,j)>a(j)|rbmxc(k,j)>a(j) if SX(j)>SY(j) && SXY(j)<=0 beta(j)=thbkfun(j); elseif SX(j)>SY(j) && SXY(j)>=0 beta(j)=thbkfun(j); elseif SY(j)>SX(j) && SXY(j)<=0 beta(j)=pi/2+thbkfun(j); elseif SY(j)>SX(j) && SXY(j)>=0 beta(j)=pi/2-thbkfun(j); end if rbmxc(k,j)>rbmx1(k,j) rbm0(k,j)=rbmxc(k,j); rbm1(k,j)=rbmx1(k,j); else rbm0(k,j)=rbmx1(k,j); rbm1(k,j)=rbmxc(k,j); end alp(k,j)=real(atanh(rbm1(k,j)./rbm0(k,j))); %Collapse Pressure Function 1: Type 1 Breakout if (STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j)) pwcffune1(k,j)=(2*cos(2*thbke+pi)*(P0(j)-C0(j).*cot(Phi(j)))-2*cosh(2*alp(k,j)).*(P0(j)-C0(j).*cot(Phi(j)))+2*SXY(j).*(-1+csc(Phi(j))).*sin(2*beta(j))... +(SXv(j)-SYv(j)).*cos(2*beta(j)).*(-1+cos(Phi(j)).*cot(Phi(j))+sin(Phi(j)))-exp(2*alp(k,j)).*(-1+cos(Phi(j)).*cot(Phi(j))+sin(Phi(j))).*((SXv(j)-SYv(j)).*cos(2*(beta(j)-thbke)-pi)+2*SXY(j).*sin(2*(beta(j)-thbke)-pi))... +(2*P0(j)-SXv(j)-SYv(j)).*(-1+csc(Phi(j))).*sinh(2*alp(k,j)))./(2*(cos(2*thbke+pi)-cosh(2*alp(k,j)))); dpwcffune1(k,j)=diff(pwcffune1(k,j),thbke); %First derivative of the collapse pressure function 1 dpwc1test=@(thbke)eval(dpwcffune1(k,j)); options1=optimoptions('fsolve','Display','iter'); %thbkfune1(j)=real(fsolve(dpwc1test,1,options1)); %thbkfune1(j)=double(vpasolve(dpwcffune1(j)==0,thbke)); thbkfune1(k,j)=0; elseif (SZwbk(k,j)>STwbk(k,j) && STwbk(k,j)>SRwbk(k,j)) %Collapse Pressure Function 2: Type 3 Breakout pwcffune2(k,j)=(exp(alp(k,j)).*(-P0(j)-SZZ(j)+SXv(j).*v(j)+SYv(j).*v(j)-exp(4*alp(k,j)).*(P0(j)+SZZ(j)-4*P0(j).*v(j)+v(j).*(SXv(j)+SYv(j)))... +2*C0(j).*cos(Phi(j))-(P0(j)-SZZ(j)+(SXv(j)+SYv(j)).*v(j)).*sin(Phi(j))... -2*exp(2*alp(k,j)).*((SXv(j)-SYv(j)).*v(j).*cos(2*beta(j)).*(-1+sin(Phi(j)))+2*SXY(j).*v(j).*sin(2*beta(j)).*(-1+sin(Phi(j)))... -cos(2*thbkn+pi)*(P0(j)+SZZ(j)-2*P0(j).*v(j)-2*C0(j).*cos(Phi(j))+(P0(j)-SZZ(j)+2*P0(j).*v(j)).*sin(Phi(j))))... +exp(4*alp(k,j)).*(2*C0(j).*cos(Phi(j))+2*(SXv(j)-SYv(j)).*v(j).*cos(2*(beta(j)-thbkn)-pi).*(-1+sin(Phi(j)))... +(-P0(j)+SZZ(j)+(-4*P0(j)+SXv(j)+SYv(j)).*v(j)).*sin(Phi(j))+4*SXY(j).*v(j).*(-1+sin(Phi(j))).*sin(2*(beta(j)-thbkn)-pi))))./(2*(cos(2*thbkn+pi)-cosh(2*alp(k,j))).*(1-2*v(j)+sin(Phi(j))+2*v(j).*sin(Phi(j)))); dpwcffune2(k,j)=diff(pwcffune2(k,j),thbkn); %First derivative of the collapse pressure function 2 dpwc2test=@(thbkn)eval(dpwcffune2(k,j)); options2=optimoptions('fsolve','Display','iter'); thbkfune1(k,j)=real(fsolve(dpwc2test,1,options2)); end else if SX(j)>SY(j) && SXY(j)<=0 beta(j)=thbkfun(j); elseif SX(j)>SY(j) && SXY(j)>=0 beta(j)=thbkfun(j); elseif SY(j)>SX(j) && SXY(j)<=0 beta(j)=pi/2+thbkfun(j); elseif SY(j)>SX(j) && SXY(j)>=0 beta(j)=pi/2-thbkfun(j); end if rbmxc(k,j)>rbmx1(k,j) rbm0(k,j)=rbmxc(k,j); rbm1(k,j)=rbmx1(k,j); else rbm0(k,j)=rbmx1(k,j); rbm1(k,j)=rbmxc(k,j); end alp(k,j)=pi/4; %Collapse Pressure Function 1: Type 1 Breakout if (STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j)) pwcffune1(k,j)=(2*cos(2*thbke+pi)*(P0(j)-C0(j).*cot(Phi(j)))-2*cosh(2*alp(k,j)).*(P0(j)-C0(j).*cot(Phi(j)))+2*SXY(j).*(-1+csc(Phi(j))).*sin(2*beta(j))... +(SXv(j)-SYv(j)).*cos(2*beta(j)).*(-1+cos(Phi(j)).*cot(Phi(j))+sin(Phi(j)))-exp(2*alp(k,j)).*(-1+cos(Phi(j)).*cot(Phi(j))+sin(Phi(j))).*((SXv(j)-SYv(j)).*cos(2*(beta(j)-thbke)-pi)+2*SXY(j).*sin(2*(beta(j)-thbke)-pi))... +(2*P0(j)-SXv(j)-SYv(j)).*(-1+csc(Phi(j))).*sinh(2*alp(k,j)))./(2*(cos(2*thbke+pi)-cosh(2*alp(k,j)))); dpwcffune1(k,j)=diff(pwcffune1(k,j),thbke); %First derivative of the collapse pressure function 1 dpwc1test=@(thbke)eval(dpwcffune1(k,j)); options1=optimoptions('fsolve','Display','iter'); %thbkfune1(j)=real(fsolve(dpwc1test,1,options1)); %thbkfune1(j)=double(vpasolve(dpwcffune1(j)==0,thbke)); thbkfune1(k,j)=0; elseif (SZwbk(k,j)>STwbk(k,j) && STwbk(k,j)>SRwbk(k,j)) %Collapse Pressure Function 2: Type 3 Breakout pwcffune2(k,j)=(exp(alp(k,j)).*(-P0(j)-SZZ(j)+SXv(j).*v(j)+SYv(j).*v(j)-exp(4*alp(k,j)).*(P0(j)+SZZ(j)-4*P0(j).*v(j)+v(j).*(SXv(j)+SYv(j)))... +2*C0(j).*cos(Phi(j))-(P0(j)-SZZ(j)+(SXv(j)+SYv(j)).*v(j)).*sin(Phi(j))... -2*exp(2*alp(k,j)).*((SXv(j)-SYv(j)).*v(j).*cos(2*beta(j)).*(-1+sin(Phi(j)))+2*SXY(j).*v(j).*sin(2*beta(j)).*(-1+sin(Phi(j)))... -cos(2*thbkn+pi)*(P0(j)+SZZ(j)-2*P0(j).*v(j)-2*C0(j).*cos(Phi(j))+(P0(j)-SZZ(j)+2*P0(j).*v(j)).*sin(Phi(j))))... +exp(4*alp(k,j)).*(2*C0(j).*cos(Phi(j))+2*(SXv(j)-SYv(j)).*v(j).*cos(2*(beta(j)-thbkn)-pi).*(-1+sin(Phi(j)))... +(-P0(j)+SZZ(j)+(-4*P0(j)+SXv(j)+SYv(j)).*v(j)).*sin(Phi(j))+4*SXY(j).*v(j).*(-1+sin(Phi(j))).*sin(2*(beta(j)-thbkn)-pi))))./(2*(cos(2*thbkn+pi)-cosh(2*alp(k,j))).*(1-2*v(j)+sin(Phi(j))+2*v(j).*sin(Phi(j)))); dpwcffune2(k,j)=diff(pwcffune2(k,j),thbkn); %First derivative of the collapse pressure function 2 dpwc2test=@(thbkn)eval(dpwcffune2(k,j)); options2=optimoptions('fsolve','Display','iter'); thbkfune1(k,j)=real(fsolve(dpwc2test,1,options2)); end end %Breakout Angle if SX(j)>SY(j) && SXY(j)<=0 thbkso(k,j)=pi/2-thbkfune1(k,j); elseif SX(j)>SY(j) && SXY(j)>=0 thbkso(k,j)=pi/2+thbkfune1(k,j); elseif SY(j)>SX(j) && SXY(j)<=0 thbkso(k,j)=pi-thbkfune1(k,j); elseif SY(j)>SX(j) && SXY(j)>=0 thbkso(k,j)=thbkfune1(k,j); end %Stresses around Wellbore wall at Maximum Compression Direction Seew1(k,j)=pwmid(k,j)-P0(j); Snnw1(k,j)=(pwmid(k,j)-P0(j))+(SYv(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*cos(2*beta(j))-cos(2*(beta(j)-thbkso(k,j))-pi))+sinh(2*alp(k,j)))... +SXv(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*cos(2*beta(j))-cos(2*(beta(j)-thbkso(k,j))-pi))+sinh(2*alp(k,j)))... +SXY(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*sin(2*beta(j))-sin(2*(beta(j)-thbkso(k,j))-pi))+sinh(2*alp(k,j)))... +SXY(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*sin(2*beta(j))-sin(2*(beta(j)-thbkso(k,j))-pi))+sinh(2*alp(k,j))))./(cosh(2*alp(k,j))-cos(2*thbkso(k,j))); SZw1(k,j)=SZZ(j)+v(j).*(Seew1(k,j)+Snnw1(k,j)); %Stresses around Wellbore wall at minimum compression direction thbksom(k,j)=pi/2+thbkso(k,j); Seew2(k,j)=pwmid(k,j)-P0(j); Snnw2(k,j)=(pwmid(k,j)-P0(j))+(SYv(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*cos(2*beta(j))-cos(2*(beta(j)-thbksom(k,j))-pi))+sinh(2*alp(k,j)))... +SXv(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*cos(2*beta(j))-cos(2*(beta(j)-thbksom(k,j))-pi))+sinh(2*alp(k,j)))... +SXY(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*sin(2*beta(j))-sin(2*(beta(j)-thbksom(k,j))-pi))+sinh(2*alp(k,j)))... +SXY(j).*(exp(2*alp(k,j)).*(exp(-alp(k,j)).*sin(2*beta(j))-sin(2*(beta(j)-thbksom(k,j))-pi))+sinh(2*alp(k,j))))./(cosh(2*alp(k,j))-cos(2*thbksom(k,j))); SZw2(k,j)=SZZ(j)+v(j).*(Seew2(k,j)+Snnw2(k,j)); if (Snnw1(k,j)>SZw1(k,j) && SZw1(k,j)>Seew1(k,j)) %Type 1 Breakout pwc(k,j)=(2*cos(2*thbkso(k,j))*(P0(j)-C0(j).*cot(Phi(j)))-2*cosh(2*alp(k,j)).*(P0(j)-C0(j).*cot(Phi(j)))+2*SXY(j).*(-1+csc(Phi(j))).*sin(2*beta(j))... +(SXv(j)-SYv(j)).*cos(2*beta(j)).*(-1+cos(Phi(j)).*cot(Phi(j))+sin(Phi(j)))-exp(2*alp(k,j)).*(-1+cos(Phi(j)).*cot(Phi(j))+sin(Phi(j))).*((SXv(j)-SYv(j)).*cos(2*(beta(j)-thbkso(k,j)))+2*SXY(j).*sin(2*(beta(j)-thbkso(k,j))))... +(2*P0(j)-SXv(j)-SYv(j)).*(-1+csc(Phi(j))).*sinh(2*alp(k,j)))./(2*(cos(2*thbkso(k,j))-cosh(2*alp(k,j)))); elseif (SZw1(k,j)>Snnw1(k,j) && Snnw1(k,j)>Seew1(k,j)) %Type 3 Breakout pwc(k,j)=((1-v(j)).*sinh(2*alp(k,j)).*(2*P0(j)-SXv(j)-SYv(j))+2*v(j).*(P0(j)*cos(2*thbkso(k,j))-SXY(j).*sin(2*beta(j)))-(P0(j)+SZZ(j))*cos(2*thbkso(k,j))+cosh(2*alp(k,j)).*(-2*v(j).*P0(j)+P0(j)+SZZ(j))... +exp(2*alp(k,j)).*(v(j)-1).*((SXv(j)-SYv(j)).*cos(2*(beta(j)-thbkso(k,j)))+2*SXY(j).*sin(2*(beta(j)-thbkso(k,j))))+(1-v(j)).*cos(2*beta(j)).*(SXv(j)+SYv(j))+2*SXY(j).*sin(2*beta(j)))./((1-2*v(j)).*(cosh(2*alp(k,j))-cos(2*thbkso(k,j)))); end end end % Stresses around the elliptical wellbore in Elliptical Coordinates e=sym('e','real'); n=sym('n','real'); ph=(e+n*1i); for k=1:Tt for j=nst:nend %Stresses due to SYv Loading only equation48(k,j)=0.25*(SYv(j).*(exp(2*alp(k,j)).*cos(2*beta(j))+(1-exp(2*alp(k,j)+2i*beta(j)))*coth(ph))... +SYv(j).*(exp(2*alp(k,j)).*cos(2*beta(j))+(1-exp(2*alp(k,j)-2i*beta(j))).*coth(2*alp(k,j)-(ph)))... +SYv(j).*coth(2*alp(k,j)-(ph)).*((1-exp(2*alp(k,j)+2*1i*beta(j))).*(csch(ph))^2)... +SYv(j).*csch(2*alp(k,j)-(ph)).*(-(cosh(2*alp(k,j))-cos(2*beta(j)))*coth(ph)*csch(ph)+2*exp(2*alp(k,j)).*cosh(2*((ph)-alp(k,j)-beta(j)*1i))*csch(ph)... -exp(2*alp(k,j)).*sinh(2*((ph)-alp(k,j)-beta(j)*1i))*coth(ph)*csch(ph))); %In-plane shear stress due to SYv loading sen1(k,j)=-imag(expand(equation48(k,j))); %Hyperbolic stress due to SYv see1(k,j)=real(expand(equation48(k,j))); %Hoops or elliptical stress due to SYv equation9a(k,j)=SYv(j).*(exp(2*alp(k,j)).*cos(2*beta(j))+(1-exp(2*alp(k,j)+2i*beta(j)))*coth(ph)); snn1(k,j)=real(expand(equation9a(k,j)))-see1(k,j); %Stresses due to SXv Loading Only equation52(k,j)=0.25*(SXv(j).*(exp(2*alp(k,j)).*cos(2*beta(j)+pi)+(1-exp(2*alp(k,j)+2*1i*(beta(j)+pi/2))).*coth(ph))... +SXv(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+pi/2))+(1-exp(2*alp(k,j)-2*1i*(beta(j)+pi/2))).*coth(2*alp(k,j)-ph))... +SXv(j).*coth(2*alp(k,j)-ph).*((1-exp(2*alp(k,j)+2*1i*(beta(j)+pi/2))).*(csch(ph))^2)... +SXv(j).*csch(2*alp(k,j)-ph).*(-(cosh(2*alp(k,j))-cos(2*(beta(j)+pi/2)))*coth(ph)*csch(ph)... +2*exp(2*alp(k,j)).*cosh(2*(ph-alp(k,j)-(beta(j)+pi/2)*1i))*csch(ph)... -exp(2*alp(k,j)).*sinh(2*(ph-alp(k,j)-(beta(j)+pi/2).*1i))*coth(ph)*csch(ph))); %In-plane shear stress due to SXv loading sen2(k,j)=-imag(expand(equation52(k,j))); %Hyperbolic stress due to SXv loading see2(k,j)=real(expand(equation52(k,j))); %Hoops stress due to SXv loading equation49(k,j)=SXv(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+pi/2))+(1-exp(2*alp(k,j)+2i*(beta(j)+pi/2)))*coth(ph)); snn2(k,j)=real(expand(equation49(k,j)))-see2(k,j); %Stresses due to SXY+ Loading equation64a(k,j)=0.25*(SXY(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+pi/4))+(1-exp(2*alp(k,j)+2*1i*(beta(j)+pi/4))).*coth(ph))... +SXY(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+pi/4))+(1-exp(2*alp(k,j)-2*1i*(beta(j)+pi/4))).*coth(2*alp(k,j)-ph))... +SXY(j).*coth(2*alp(k,j)-ph).*((1-exp(2*alp(k,j)+2*1i*(beta(j)+pi/4))).*(csch(ph))^2)... +SXY(j).*csch(2*alp(k,j)-ph).*(-(cosh(2*alp(k,j))-cos(2*(beta(j)+pi/4)))*coth(ph)*csch(ph)... +2*exp(2*alp(k,j)).*cosh(2*(ph-alp(k,j)-(beta(j)+pi/4)*1i))*csch(ph)... -exp(2*alp(k,j)).*sinh(2*(ph-alp(k,j)-(beta(j)+pi/4).*1i))*coth(ph)*csch(ph))); %In-plane shear stress due to SXY+ loading sen3(k,j)=-imag(expand(equation64a(k,j))); %Hyperbolic stress due to SXY+ loading see3(k,j)=real(expand(equation64a(k,j))); %Hoops stress due to SXY+ loading equation64ab(k,j)=SXY(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+pi/4))+(1-exp(2*alp(k,j)+2i*(beta(j)+pi/4)))*coth(ph)); snn3(k,j)=real(expand(equation64ab(k,j)))-see3(k,j); %Stresses due to SXY- Loading equation64b(k,j)=0.25*(-SXY(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+3*pi/4))+(1-exp(2*alp(k,j)+2*1i*(beta(j)+3*pi/4))).*coth(ph))... -SXY(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+3*pi/4))+(1-exp(2*alp(k,j)-2*1i*(beta(j)+3*pi/4))).*coth(2*alp(k,j)-ph))... -SXY(j).*coth(2*alp(k,j)-ph).*((1-exp(2*alp(k,j)+2*1i*(beta(j)+3*pi/4))).*(csch(ph))^2)... -SXY(j).*csch(2*alp(k,j)-ph).*(-(cosh(2*alp(k,j))-cos(2*(beta(j)+3*pi/4)))*coth(ph)*csch(ph)... +2*exp(2*alp(k,j)).*cosh(2*(ph-alp(k,j)-(beta(j)+3*pi/4)*1i))*csch(ph)... -exp(2*alp(k,j)).*sinh(2*(ph-alp(k,j)-(beta(j)+3*pi/4).*1i))*coth(ph)*csch(ph))); %In-plane shear stress due to SXY- loading sen4(k,j)=-imag(expand(equation64b(k,j))); %Hyperbolic stress due to SXY- loading see4(k,j)=real(expand(equation64b(k,j))); %Hoops stress due to SXY- loading equation64abc(k,j)=-SXY(j).*(exp(2*alp(k,j)).*cos(2*(beta(j)+3*pi/4))+(1-exp(2*alp(k,j)+2i*(beta(j)+3*pi/4)))*coth(ph)); snn4(k,j)=real(expand(equation64abc(k,j)))-see4(k,j); %Total Stresses snn(k,j)=snn1(k,j)+snn2(k,j)+snn3(k,j)+snn4(k,j)+(pwmid(k,j)-P0(j));%hoops stress see(k,j)=see1(k,j)+see2(k,j)+see3(k,j)+see4(k,j)+(pwmid(k,j)-P0(j));%hyperbolic stress sen(k,j)=sen1(k,j)+sen2(k,j)+sen3(k,j)+sen4(k,j);%in-plane shear stress SZ(k,j)=SZZ(j)+v(j).*(see(k,j)+snn(k,j)); %Vertical stress...plane strain approx. end end for k=1:Tt for j=nst:nend %Stresses in the direction of max compression for all distances from the %wellbore snnbk(k,j)=subs(snn(k,j),n,thbkso(k,j)); seebk(k,j)=subs(see(k,j),n,thbkso(k,j)); senbk(k,j)=subs(sen(k,j),n,thbkso(k,j)); szbk(k,j)=subs(SZ(k,j),n,thbkso(k,j)); %Stresses in the direction of min. compression for all distances from the %wellbore thbksom(k,j)=thbkso(k,j); snnbkm(k,j)=subs(snn(k,j),n,(thbksom(k,j))); seebkm(k,j)=subs(see(k,j),n,(thbksom(k,j))); senbkm(k,j)=subs(sen(k,j),n,(thbksom(k,j))); szbkm(k,j)=subs(SZ(k,j),n,(thbksom(k,j))); %At the wellbore wall, the stresses in the direction of maximum compression %are snnbkw(k,j)=double(subs(snnbk(k,j),e,alp(k,j))); seebkw(k,j)=double(subs(seebk(k,j),e,alp(k,j))); senbkw(k,j)=double(subs(senbk(k,j),e,alp(k,j))); SZwbk(k,j)=double(subs(szbk(k,j),e,alp(k,j))); %At the wellbore wall, the stresses in the direction of minimum compression snnbkwm(k,j)=double(subs(snnbkm(k,j),e,alp(k,j))); seebkwm(k,j)=double(subs(seebkm(k,j),e,alp(k,j))); senbkwm(k,j)=double(subs(senbkm(k,j),e,alp(k,j))); SZwbkm(k,j)=double(subs(szbkm(k,j),e,alp(k,j))); end end %Estimating Breakout Depth in the direction of maximum compression; noting %that the wellbore geometry is now assumed to be elliptical in shape for k=1:Tt for j=nst:nend if rbmx1(k,j)>a(j)|rbmxc(k,j)>a(j) if pwmid(k,j)<=pwc(k,j) if (snnbkw(k,j)>SZwbk(k,j) && SZwbk(k,j)>seebkw(k,j)) %Type 1 breakout sigma1A(k,j)=0.5*(snnbk(k,j)+seebk(k,j))+0.5*((snnbk(k,j)-seebk(k,j)).^2+4*senbk(k,j).^2).^0.5; sigma3A(k,j)=0.5*(snnbk(k,j)+seebk(k,j))-0.5*((snnbk(k,j)-seebk(k,j)).^2+4*senbk(k,j).^2).^0.5; elseif (SZwbk(k,j)>snnbkw(k,j) && snnbkw(k,j)>seebkw(k,j)) sigma1A(k,j)=szbk(k,j); sigma3A(k,j)=0.5*(snnbk(k,j)+seebk(k,j))-0.5*((snnbk(k,j)-seebk(k,j)).^2+4*senbk(k,j).^2).^0.5; end C0fun(k,j)=0.5*((sigma1A(k,j).*((1+mu(j).^2).^0.5-mu(j)))-(sigma3A(k,j).*((1+mu(j).^2).^0.5+mu(j)))); %Cohesive Strength Function C0funtest=@(e)eval(C0fun(k,j)-C0(j)); options=optimoptions('fsolve','Display','iter'); e1(k,j)=real(fsolve(C0funtest,0.7,options)); %esoln{j}=(double(vpasolve(C0fun(j)==C0(j),e)));% using matlab solver %e1(j)=esoln{j}(1); if e1(k,j)<0 rbmxe(k,j)=a(j)./rD1(k,j); else f1(k,j)=sqrt((a(j)./rD1(k,j)).^2-(a(j)./rD1m(k,j)).^2); %focal length of the confocal ellipses rbmxp(k,j)=f1(k,j)./sqrt(1-(tanh(e1(k,j))).^2); rbmxe(k,j)= rbmxp(k,j)./cos(thbkfun(j)-thbkfune1(k,j)); %breakout depth in the direction of max. compression if rbmxe(k,j)<=0 rbmxe(k,j)=rbmx1(k,j); end end else e1(k,j)=alp(k,j); rbmxe(k,j)=rbmx1(k,j); end else e1(k,j)=alp(k,j); rbmxe(k,j)=rbmx1(k,j); end end end %Estimating Breakout Depth in the direction of minimum compression; for k=1:Tt for j=nst:nend if rbmx1(k,j)>a(j)|rbmxc(k,j)>a(j) if pwmid(k,j)<=pwc(k,j) if (snnbkwm(k,j)>SZwbkm(k,j) && SZwbkm(k,j)>seebkwm(k,j)) %Type 1 breakout sigma1Am(k,j)=0.5*(snnbkm(k,j)+seebkm(k,j))+0.5*((snnbkm(k,j)-seebkm(k,j)).^2+4*senbkm(k,j).^2).^0.5; sigma3Am(k,j)=0.5*(snnbkm(k,j)+seebkm(k,j))-0.5*((snnbkm(k,j)-seebkm(k,j)).^2+4*senbkm(k,j).^2).^0.5; %sigma3Amr(j)=real(expand(sigma3Am(j))); elseif (SZwbkm(k,j)>snnbkwm(k,j) && snnbkwm(k,j)>seebkwm(k,j))%Type 3 breakout sigma1Am(k,j)=szbkm(k,j); sigma3Am(k,j)=0.5*(snnbkm(k,j)+seebkm(k,j))-0.5*((snnbkm(k,j)-seebkm(k,j)).^2+4*senbkm(k,j).^2).^0.5; %sigma3Amr(j)=real(expand(sigma3Am(j))); end C0funm(k,j)=0.5*((sigma1Am(k,j).*((1+mu(j).^2).^0.5-mu(j)))-(sigma3Am(k,j).*((1+mu(j).^2).^0.5+mu(j)))); %Cohesive Strength Function %C0funmr(j)=real(expand(C0funm(j))); C0funmtest=@(e)eval(C0funm(k,j)-C0(j)); options=optimoptions('fsolve','Display','iter'); e1m(k,j)=real(fsolve(C0funmtest,0.7,options)); %esolnm{j}=(double(vpasolve(C0funm(j)==C0(j),e)));% using matlab solver %e1m(j)=esolnm{j}(1); if e1m(k,j)<0 rbmxce(k,j)=a(j)./rD1m(k,j); else f1(k,j)=sqrt((a(j)./rD1(k,j)).^2-(a(j)./rD1m(k,j)).^2); %focal length of the confocal ellipses rbmxpc(k,j)=f1(k,j)./sqrt(1-(tanh(e1m(k,j))).^2); rbmxce(k,j)= rbmxpc(k,j)./cos(thbkfun(j)-thbkfune1(k,j)); %breakout depth in the direction of max. compression if rbmxce(k,j)<=0 rbmxce(k,j)=rbmxc(k,j); end end else e1m(k,j)=alp(k,j); rbmxce(k,j)=0; end else e1m(k,j)=alp(k,j); rbmxce(k,j)=0; end end end %Breakout Width in the direction of maximum compression...elliptic wellbore for k=1:Tt for j=nst:nend if rbmx1(k,j)>a(j)|rbmxc(k,j)>a(j) if pwmid(k,j)<=pwc(j) %Critical Hoops Stress STTc(k,j)=(2*C0(j).*cos(Phi(j)))./(1-sin(Phi(j)))+((pwc(k,j)-P0(j)).*(1+sin(Phi(j))))./(1-sin(Phi(j))); if (snnbkw(k,j)>SZwbk(k,j) && SZwbk(k,j)>seebkw(k,j)) %Type 1 breakout D1(k,j)=STTc(k,j)-(pwc(k,j)-P0(j))+exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*cos(2*beta(j))+2*exp(2*alp(k,j)).*SXY(j).*sin(2*beta(j)); B1(k,j)=exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*sin(2*beta(j))-2*exp(2*alp(k,j)).*SXY(j).*cos(2*beta(j)); R1(k,j)=sqrt(B1(k,j).^2+D1(k,j).^2); omega1(k,j)=atan(B1(k,j)./D1(k,j)); C1(k,j)=STTc(k,j).*cosh(2*alp(k,j))+(SXv(j)-SYv(j)).*cos(2*beta(j))-(pwc(k,j)-P0(j)).*cosh(2*alp(k,j))+2*SXY(j).*sin(2*beta(j))+(2*P0(j)-SXv(j)-SYv(j)).*sinh(2*alp(k,j)); thetabkte1(k,j)=double(real(0.5*(acos(C1(k,j)./R1(k,j))-omega1(k,j)))); if thetabkte1(k,j)<0 thetabkte(k,j)=pi+thetabkte1(k,j); elseif thetabkte1(k,j)>2*pi thetabkte(k,j)=2*pi; else thetabkte(k,j)=thetabkte1(k,j); end elseif (SZwbk(k,j)>snnbkw(k,j) && snnbkw(k,j)>seebkw(k,j))%Type 3 breakout D1(k,j)=STTc(k,j)-(v(j)+1).*(pwc(k,j)-P0(j))-SZZ(j)+v(j).*exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*cos(2*beta(j))+2*v(j).*exp(2*alp(k,j)).*SXY(j).*sin(2*beta(j)); B1(k,j)=v(j).*exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*sin(2*beta(j))-2*v(j).*exp(2*alp(k,j)).*SXY(j).*cos(2*beta(j)); R1(k,j)=sqrt(B1(k,j).^2+D1(k,j).^2); omega1(k,j)=atan(B1(k,j)./D1(k,j)); C1(k,j)=STTc(k,j).*cosh(2*alp(k,j))-SZZ(j).*cosh(2*alp(k,j))+v(j).*(SXv(j)-SYv(j)).*cos(2*beta(j))-v(j).*(pwc(k,j)-P0(j)).*cosh(2*alp(k,j))+2*v(j).*SXY(j).*sin(2*beta(j))+v(j).*(2*P0(j)-SXv(j)-SYv(j)).*sinh(2*alp(k,j)); thetabkte1(k,j)=double(real(0.5*(acos(C1(k,j)./R1(k,j))-omega1(k,j)))); if thetabkte1(k,j)<0 thetabkte(k,j)=pi+thetabkte1(k,j); elseif thetabkte1(k,j)>2*pi thetabkte(k,j)=2*pi; else thetabkte(k,j)=thetabkte1(k,j); end end end end end end %Breakout Width in the direction of minimum compression...elliptic wellbore for k=1:Tt for j=nst:nend if rbmx1(k,j)>a(j)|rbmxc(k,j)>a(j) if pwmid(k,j)<=pwc(k,j) pwcme(k,j)=subs(pwc(k,j),n,(thbksom(k,j))); %Critical Hoops Stress STTcm(k,j)=(2*C0(j).*cos(Phi(j)))./(1-sin(Phi(j)))+((pwcme(k,j)-P0(j)).*(1+sin(Phi(j))))./(1-sin(Phi(j))); if (snnbkwm(k,j)>SZwbkm(k,j) && SZwbkm(k,j)>seebkwm(k,j)) %Type 1 breakout D1m(k,j)=STTc(k,j)-(pwmid(k,j)-P0(j))+exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*cos(2*beta(j))+2*exp(2*alp(k,j)).*SXY(j).*sin(2*beta(j)); B1m(k,j)=exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*sin(2*beta(j))-2*exp(2*alp(k,j)).*SXY(j).*cos(2*beta(j)); R1m(k,j)=sqrt(B1m(k,j).^2+D1m(k,j).^2); omega1m(k,j)=atan(B1m(k,j)./D1m(k,j)); C1m(k,j)=STTcm(k,j).*cosh(2*alp(k,j))+(SXv(j)-SYv(j)).*cos(2*beta(j))-(pwmid(k,j)-P0(j)).*cosh(2*alp(k,j))+2*SXY(j).*sin(2*beta(j))+(2*P0(j)-SXv(j)-SYv(j)).*sinh(2*alp(k,j)); thetabkte1m(k,j)=double(real(0.5*(acos(C1m(k,j)./R1m(k,j))-omega1m(k,j)))); if thetabkte1m(k,j)<0 thetabktem(k,j)=2*pi+thetabkte1m(k,j); elseif thetabkte1m(k,j)>2*pi thetabktem(k,j)=2*pi; else thetabktem(k,j)=thetabkte1m(k,j); end elseif (SZwbkm(k,j)>snnbkwm(k,j) && snnbkwm(k,j)>seebkwm(k,j))%Type 3 breakout D1m(k,j)=STTcm(k,j)-(v(j)+1).*(pwmid(k,j)-P0(j))-SZZ(j)+v(j).*exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*cos(2*beta(j))+2*v(j).*exp(2*alp(k,j)).*SXY(j).*sin(2*beta(j)); B1m(k,j)=v(j).*exp(2*alp(k,j)).*(SXv(j)-SYv(j)).*sin(2*beta(j))-2*v(j).*exp(2*alp(k,j)).*SXY(j).*cos(2*beta(j)); R1m(k,j)=sqrt(B1m(k,j).^2+D1m(k,j).^2); omega1m(k,j)=atan(B1m(k,j)./D1m(k,j)); C1m(k,j)=STTcm(k,j).*cosh(2*alp(k,j))-SZZ(j).*cosh(2*alp(k,j))+v(j).*(SXv(j)-SYv(j)).*cos(2*beta(j))-v(j).*(pwmid(k,j)-P0(j)).*cosh(2*alp(k,j))+2*v(j).*SXY(j).*sin(2*beta(j))+v(j).*(2*P0(j)-SXv(j)-SYv(j)).*sinh(2*alp(k,j)); thetabkte1m(k,j)=double(real(0.5*(acos(C1m(k,j)./R1m(k,j))-omega1m(k,j)))); if thetabkte1m(k,j)<0 thetabktem(k,j)=2*pi+thetabkte1m(k,j); elseif thetabkte1m(k,j)>2*pi thetabktem(k,j)=2*pi; else thetabktem(k,j)=thetabkte1m(k,j); end end end end end end % Equivalent Radius of the Wellbore and Pseudo Breakout Width for k=1:Tt %time for j=nst:nend %position if rbmx1(k,j)>a(j)|rbmxc(k,j)>a(j) %first episode from circular wellbore rbktm(j)=max(rbmxc(k,j)); rbmx(j)=max(rbmx1(k,j)); elseif rbmxe(k,j)>rbmx1(k,j)|rbmxce(k,j)>rbmxc(k,j)%breakout from elliptical wellbore rbktm(j)=max(rbmxce(k,j)); rbmx(j)=max(rbmxe(k,j)); end rwr(j)=sqrt(rbmx(j).*rbktm(j));%Using Equivalent Area Hypothesis ...in inches rbkt_art(j) = 0.0254*(rwr(j)-a(j)+rbkt_int(j)); % Pseudo Differential Breakout Depth...Comsol Link... in meters end end % %Defining the Collapse Pressure for Use in Comsol Deskstop...unit in Psi % pwc1=pwc(1,1);pwc2=pwc(2,1);pwc3=pwc(3,1);pwc4=pwc(4,1);pwc5=pwc(5,1);pwc6=pwc(6,1);pwc7=pwc(7,1);pwc8=pwc(8,1);pwc9=pwc(9,1);pwc10=pwc(10,1); % pwc11=pwc(11,1);pwc12=pwc(12,1);pwc13=pwc(13,1);pwc14=pwc(14,1);pwc15=pwc(15,1);pwc16=pwc(16,1);pwc17=pwc(17,1);pwc18=pwc(18,1);pwc19=pwc(19,1);pwc20=pwc(20,1); % pwc21=pwc(21,1);pwc22=pwc(22,1);pwc23=pwc(23,1);pwc24=pwc(24,1);pwc25=pwc(25,1);pwc26=pwc(26,1);pwc27=pwc(27,1);pwc28=pwc(28,1); % %Linking the Critical Collapse Pressure in Matlab with Comsol in Psi % model.param.set('pwc1',pwc1);model.param.set('pwc2',pwc2);model.param.set('pwc3',pwc3);model.param.set('pwc4',pwc4); % model.param.set('pwc5',pwc5);model.param.set('pwc6',pwc6);model.param.set('pwc7',pwc7);model.param.set('pwc8',pwc8); % model.param.set('pwc9',pwc9);model.param.set('pwc10',pwc10);model.param.set('pwc11',pwc11);model.param.set('pwc12',pwc12); % model.param.set('pwc13',pwc13);model.param.set('pwc14',pwc14);model.param.set('pwc15',pwc15);model.param.set('pwc16',pwc16); % model.param.set('pwc17',pwc17);model.param.set('pwc18',pwc18);model.param.set('pwc19',pwc19);model.param.set('pwc20',pwc20); % model.param.set('pwc21',pwc21);model.param.set('pwc22',pwc22);model.param.set('pwc23',pwc23);model.param.set('pwc24',pwc24); % model.param.set('pwc25',pwc25);model.param.set('pwc26',pwc26);model.param.set('pwc27',pwc27);model.param.set('pwc28',pwc28); % %Linking the Mid Average Pressure Values to COMSOL in Psi % model.param.set('pwmid1',pwmid1);model.param.set('pwmid2',pwmid2);model.param.set('pwmid3',pwmid3);model.param.set('pwmid4',pwmid4); % model.param.set('pwmid5',pwmid5);model.param.set('pwmid6',pwmid6);model.param.set('pwmid7',pwmid7);model.param.set('pwmid8',pwmid8); % model.param.set('pwmid9',pwmid9);model.param.set('pwmid10',pwmid10);model.param.set('pwmid11',pwc11);model.param.set('pwmid12',pwmid12); % model.param.set('pwmid13',pwmid13);model.param.set('pwmid14',pwmid14);model.param.set('pwmid15',pwmid15);model.param.set('pwmid16',pwmid16); % model.param.set('pwmid17',pwmid17);model.param.set('pwmid18',pwmid18);model.param.set('pwc19',pwc19);model.param.set('pwc20',pwc20); % model.param.set('pwc21',pwc21);model.param.set('pwc22',pwc22);model.param.set('pwc23',pwc23);model.param.set('pwc24',pwc24); % model.param.set('pwc25',pwc25);model.param.set('pwc26',pwc26);model.param.set('pwc27',pwc27);model.param.set('pwc28',pwc28); %Updating Breakout Depth in Comsol in meters rbkt_art1=rbkt_art(1);rbkt_art2=rbkt_art(2);rbkt_art3=rbkt_art(3);rbkt_art4=rbkt_art(4);rbkt_art5=rbkt_art(5);rbkt_art6=rbkt_art(6);rbkt_art7=rbkt_art(7); rbkt_art8=rbkt_art(8);rbkt_art9=rbkt_art(9);rbkt_art10=rbkt_art(10); rbkt_art11=rbkt_art(11);rbkt_art12=rbkt_art(12);rbkt_art13=rbkt_art(13);rbkt_art14=rbkt_art(14);rbkt_art15=rbkt_art(15); rbkt_art16=rbkt_art(16);rbkt_art17=rbkt_art(17);rbkt_art18=rbkt_art(18); % rbkt_art19=rbkt_art(19);rbkt_art20=rbkt_art(20);rbkt_art21=rbkt_art(21); % % rbkt_art22=rbkt_art(22);rbkt_art23=rbkt_art(23);rbkt_art24=rbkt_art(24);rbkt_art25=rbkt_art(25);rbkt_art26=rbkt_art(26);rbkt_art27=rbkt_art(27);rbkt_art28=rbkt_art(28); % model.param.set('tL1',rbkt_art1);model.param.set('tL2',rbkt_art2);model.param.set('tL3',rbkt_art3);model.param.set('tL4',rbkt_art4);model.param.set('tL5',rbkt_art5); model.param.set('tL6',rbkt_art6);model.param.set('tL7',rbkt_art7);model.param.set('tL8',rbkt_art8);model.param.set('tL9',rbkt_art9);model.param.set('tL10',rbkt_art10); model.param.set('tL11',rbkt_art11);model.param.set('tL12',rbkt_art12);model.param.set('tL13',rbkt_art13);model.param.set('tL14',rbkt_art14);model.param.set('tL15',rbkt_art15); model.param.set('tL16',rbkt_art16);model.param.set('tL17',rbkt_art17);model.param.set('tL18',rbkt_art18); % model.param.set('tL19',rbkt_art19);model.param.set('tL20',rbkt_art20); % model.param.set('tL21',rbkt_art21);model.param.set('tL22',rbkt_art22);model.param.set('tL23',rbkt_art23);model.param.set('tL24',rbkt_art24);model.param.set('tL25',rbkt_art25); % model.param.set('tL26',rbkt_art26);model.param.set('tL27',rbkt_art27);model.param.set('tL28',rbkt_art28); % %Post Processing...Wellbore Geometry at Producing Interval in Inches rwp1=[rbmx(1) rbmx(2) rbmx(3) rbmx(4) rbmx(5) rbmx(6) rbmx(7) rbmx(8) rbmx(9) rbx(10) rbmx(11) rbmx(12) rbmx(13) rbmx(14) rbmx(15) rbmx(16)... rbmx(17) rbmx(18)];% rbmx(19) rbmx(20) rbmx(21) rbmx(22) rbmx(23) rbmx(24) rbmx(25) rbmx(26) rbmx(27) rbmx(28)]; %half major hole axis of the openhole formations rwp2 = [rbktm(1) rbktm(2) rbktm(3) rbktm(4) rbktm(5) rbktm(6) rbktm(7) rbktm(8) rbktm(9) rbktm(10) rbktm(11) rbktm(12) rbktm(13) rbktm(14) rbktm(15)... rbktm(16) rbktm(17) rbktm(18)];% rbktm(19) rbktm(20) rbktm(21) rbktm(22) rbktm(23) rbktm(24) rbktm(25) rbktm(26) rbktm(27) rbktm(28)];%half minor hole axis of the openhole formations % ESTIMATING THE VOLUME OF SOLID INFLUX TO THE WELLBORE DURING % COLLAPSE FROM EACH LAYER for k=1:Tt for j=nst:1:nend SRw3(k,j)=pwmid(k,j)-P0(j); STw3(k,j)=(SXv(j)+SYv(j))-2*(SXv(j)-SYv(j)).*cos(2*thbkso(k,j)+pi)-4*SXY(j)*sin(2*thbkso(k,j)+pi)-pwmid(k,j)+P0(j); SZw3(k,j)=SZZ(j)-v(j).*(2*(SXv(j)-SYv(j))*cos(2*thbkso(k,j)+pi)+4*SXY(j)*sin(2*thbkso(k,j)+pi)); STZw3(k,j)=2*(-SXZ(j)*sin(thbkso(k,j)+pi)+SYZ(j)*cos(thbkso(k,j)+pi)); frac(j)=pi/4-Phi(j)/2;% Fracture Angle in Type 3 Breakout hbmxc(k,j)=(1/12)*2*(rbmxc(k,j)-a(j))./tan(frac(j)); mc(k,j)=(1/12)*(rbmxc(k,j)-a(j))./hbmxc(k,j); hbmxce(k,j)=(1/12)*2*(rbmxce(k,j)-a(j))./tan(frac(j)); mce(k,j)=(1/12)*(rbmxce(k,j)-a(j))./hbmxce(k,j); hbmx(k,j)=(1/12)*2*(rbmx1(k,j)-a(j))./tan(frac(j)); %hbmx(j)=5/12; m(k,j)=(1/12)*(rbmx1(k,j)-a(j))./hbmx(k,j); hbmxe(k,j)=(1/12)*2*(rbmxe(k,j)-a(j))./tan(frac(j)); %hbmxe(j)=5; me(k,j)=(1/12)*(rbmxe(k,j)-a(j))./hbmxe(k,j); if (STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j))||((snnbkw(k,j)>SZwbk(k,j) && SZwbk(k,j)>seebkw(k,j))) %types 1 and 3 in min and max if SXY(j)~=0 vol_sand(k,j)=(1/144)*(1-por(j)).*(thetabkt(k,j).*(rbmx1(k,j).*a(j)-a(j).^2)+thetabktm(k,j).*(rbmxc(k,j).*a(j)-a(j).^2)... +thetabkte(k,j).*(rbmxe(k,j).*a(j)-a(j).^2)+thetabktem(k,j).*(rbmxce(k,j).*a(j)-a(j).^2)).*h(j);%volume of solid grains from each layer in cubic feet elseif SXY(j)==0 if thetabkte(k,j)>thetabkt(k,j) && rbmxe(k,j)>rbmx1(k,j) vol_sand(k,j)=(1/144)*(1-por(j)).*(thetabkte(k,j).*(rbmxe(k,j).*a(j)-a(j).^2)+thetabktem(k,j).*(rbmxce(k,j).*a(j)-a(j).^2)).*h(j); elseif thetabkt(k,j)>thetabkte(k,j) && rbmxe(k,j)>rbmx1(k,j) vol_sand(k,j)=(1/144)*(1-por(j)).*(thetabkt(k,j).*(rbmx1(k,j).*a(j)-a(j).^2)+thetabktm(k,j).*(rbmxc(k,j).*a(j)-a(j).^2)... +thetabkte(k,j).*(rbmxe(k,j).*a(j)-a(j).*rbmx1(k,j))+thetabktem(k,j).*(rbmxce(k,j).*a(j)-a(j).*rbmxc(k,j))).*h(j); else vol_sand(k,j)=(1/144)*(1-por(j)).*(thetabkt(k,j).*(rbmx1(k,j).*a(j)-a(j).^2)+thetabktm(k,j).*(rbmxc(k,j).*a(j)-a(j).^2)).*h(j); end end elseif ((STwbk(k,j)>SZwbk(k,j) && SZwbk(k,j)>SRwbk(k,j)) &&(SZwbkm(k,j)>STwbkm(k,j) && STwbkm(k,j)>SRwbkm(j)))||((snnbkw(k,j)>SZwbk(k,j) && SZwbk(k,j)>seebkw(k,j))&&(SZwbkm(k,j)>snnbkwm(k,j) && snnbkwm(k,j)>seebkwm(k,j))) % type 1 in max type 3 in min if thetabkt(k,j)<=pi if SXY(j)~=0 vol_sand(k,j)=(1/144)*((1-por(j)).*(thetabkt(k,j).*(rbmx1(k,j).*a(j)-a(j).^2)).*h(j)+((1/3)*thetabktm(k,j).*(mc(j).^2.*(hbmxc(k,j).*(12*(0.5*hbmxc(k,j)).^2-3*hbmxc(k,j).^2+(hbmxc(k,j)).^2))))... +(thetabkte(k,j).*(rbmxe(k,j).*a(j)-a(j).^2)).*h(j)+((1/3)*thetabktem(k,j).*(mce(k,j).^2.*(hbmxce(k,j).*(12*(0.5*hbmxce(k,j)).^2-3*hbmxce(k,j).^2+(hbmxce(k,j)).^2)))));%volume of solid grains from each layer in cubic feet elseif SXY(j)==0 if thetabkte(k,j)>thetabkt(k,j) && rbmxce(k,j)>rbmxc(k,j) vol_sand(k,j)=(1/144)*(1-por(j)).*((thetabkte(k,j).*(rbmxe(k,j).*a(j)-a(j).^2)).*h(j)+((1/3)*thetabktem(k,j).*(mce(k,j).^2.*(hbmxce(k,j).*(12*(0.5*hbmxce(k,j)).^2-3*hbmxce(k,j).^2+(hbmxce(k,j)).^2))))); elseif thetabkt(k,j)>thetabkte(j) && rbmxce(k,j)>rbmxc(k,j) vol_sand(k,j)=(1/144)*(1-por(j)).*(thetabkt(k,j).*(rbmx1(k,j).*a(j)-a(j).^2).*h(j)+((1/3)*thetabktm(k,j).*(mc(j).^2.*(hbmxc(k,j).*(12*(0.5*hbmxc(k,j)).^2-3*hbmxc(k,j).^2+(hbmxc(k,j)).^2))))... +(thetabkte(k,j).*(rbmxe(k,j).*a(j)-a(j).*rbmx1(k,j))).*h(j)+((1/3)*thetabktem(k,j).*(mce(k,j).^2.*(hbmxce(k,j).*(12*(0.5*hbmxce(k,j)).^2-3*hbmxce(k,j).^2+(hbmxce(k,j)).^2))))); else vol_sand(k,j)=(1/144)*(1-por(j)).*(((1/3)*thetabkt(k,j).*(m(k,j).^2.*(hbmx(k,j).*(12*(0.5*hbmx(k,j)).^2-3*hbmx(k,j).^2+(hbmx(k,j)).^2))))... +((1/3)*thetabktm(k,j).*(mc(k,j).^2.*(hbmxc(k,j).*(12*(0.5*hbmxc(k,j)).^2-3*hbmxc(k,j).^2+(hbmxc(k,j)).^2))))); end end end elseif (SZwbk(k,j)>STwbk(k,j) && STwbk(k,j)>SRwbk(k,j))||(SZwbk(k,j)>snnbkw(k,j) && snnbkw(k,j)>seebkw(k,j)) %type 3 in min and max hbmx(k,j)=2*(rbmx1(k,j)-a(j))./tan(frac(j)); if SXY(j)~=0 vol_sand(k,j)=(1/144)*(1-por(j)).*(((1/3)*thetabkt(k,j).*(m(j).^2.*(hbmx(k,j).*(12*(0.5*hbmx(k,j)).^2-3*hbmx(k,j).^2+(hbmx(k,j)).^2))))... +((1/3)*thetabktm(k,j).*(mc(k,j).^2.*(hbmxc(k,j).*(12*(0.5*hbmxc(k,j)).^2-3*hbmxc(k,j).^2+(hbmxc(k,j)).^2)))))... +(1/144)*(1-por(j)).*(((1/3)*thetabkte(k,j).*(me(k,j).^2.*(hbmxe(k,j).*(12*(0.5*hbmxe(k,j)).^2-3*hbmxe(k,j).^2+(hbmxe(k,j)).^2))))... +((1/3)*thetabktme(k,j).*(mce(k,j).^2.*(hbmxce(k,j).*(12*(0.5*hbmxce(k,j)).^2-3*hbmxce(k,j).^2+(hbmxce(k,j)).^2)))));%volume of solid grains from each layer in cubic feet elseif SXY(j)==0 if thetabkte(k,j)>thetabkt(k,j) && rbmxe(k,j)>rbmx1(k,j) vol_sand(k,j)=(1/144)*(1-por(j)).*(((1/3)*thetabkte(k,j).*(me(k,j).^2.*(hbmxe(k,j).*(12*(0.5*hbmxe(k,j)).^2-3*hbmxe(k,j).^2+(hbmxe(k,j)).^2))))... +((1/3)*thetabktme(k,j).*(mce(k,j).^2.*(hbmxce(k,j).*(12*(0.5*hbmxce(k,j)).^2-3*hbmxce(k,j).^2+(hbmxce(k,j)).^2))))); elseif thetabkt(k,j)>thetabkte(k,j) && rbmxe(k,j)>rbmx1(k,j) hbmxce2(k,j)=2*(rbmxce(k,j)-rbmxc(k,j))./tan(frac(j)); mce2(k,j)=(rbmxce(k,j)-rbmxc(k,j))./hbmxce(k,j); hbmxe2(k,j)=2*(rbmxe(k,j)-rbmx1(k,j))./tan(frac(j)); me2(k,j)=(rbmxe(k,j)-rbmx1(k,j))./hbmxe(k,j); vol_sand(k,j)=(1/144)*(1-por(j)).*(((1/3)*thetabkt(k,j).*(m(k,j).^2.*(hbmx(k,j).*(12*(0.5*hbmx(k,j)).^2-3*hbmx(k,j).^2+(hbmx(k,j)).^2))))... +((1/3)*thetabktm(k,j).*(rbmxc(k,j).*mc(j).^2.*(hbmxc(k,j).*(12*(0.5*hbmxc(k,j)).^2-3*hbmxc(k,j).^2+(hbmxc(k,j)).^2)))))... +(1/144)*(1-por(j)).*(((1/3)*thetabkte(k,j).*(me2(k,j).^2.*(hbmxe2(k,j).*(12*(0.5*hbmxe2(k,j)).^2-3*hbmxe2(k,j).^2+(hbmxe2(k,j)).^2))))... +((1/3)*thetabktme(k,j).*(mce2(k,j).^2.*(hbmxce2(k,j).*(12*(0.5*hbmxce2(k,j)).^2-3*hbmxce2(k,j).^2+(hbmxce2(k,j)).^2))))); else vol_sand(k,j)=(1/144)*(1-por(j)).*(((1/3)*thetabkt(k,j).*(m(k,j).^2.*(hbmx(k,j).*(12*(0.5*hbmx(k,j)).^2-3*hbmx(k,j).^2+(hbmx(k,j)).^2))))... +((1/3)*thetabktm(k,j).*(rbmxc(k,j).*mc(k,j).^2.*(hbmxc(k,j).*(12*(0.5*hbmxc(k,j)).^2-3*hbmxc(k,kj).^2+(hbmxc(k,j)).^2))))); end end end vol_sands(j)=max(vol_sand(k,j)); end end Vol_sand_total = sum(vol_sands); %total solid influx %ESTIMATING THE VOLUME OF THE WELLBORE-RISER SYSTEM %Volume of each open hole layers for j=nst:1:nend vol_well(j) = max(vol_sand(k,j))+(pi-min(thetabkt(k,j))).*a(j).*h(j)+(pi-min(thetabktm(k,j))).*a(j).*h(j); %in cubic feet Volwel_ophl = sum(vol_well); end for k=1:1:nc volc_well(k)=(1/144)*pi*rc(k).*rc(k).*hc(k); %non-openhole layers with different casing ID in cubic feet Volwel_csd = sum(volc_well); end Vol_Well_Total =Volwel_ophl+Volwel_csd; % total wellbore-riser volume in cubic feet %Approximated Fluid Properties rhostar=rhof*(1+Vol_sand_total/Vol_Well_Total); %in g/cm3 mustar = mufw/rhof*((1-Vol_sand_total/Vol_Well_Total)*rhostar+(Vol_sand_total/Vol_Well_Total)*rhos);% in Pa.s % Chien 1994 Empirical Model for Settling Velocity in cm/s vst = 120*(mustar/(dp*rhostar))*((1+0.0727*dp*(rhos/rhostar-1)*(dp*rhostar/mustar)^2)^0.5-1); % Solid Settling Mass Flux for j=nst:1:nend for k=1:Tt if pwmid(k,j)<=pwc(k,j) area_sand(j)=(6.4516)*max(thetabkt(k,j)).*(max(rbmx(k,j)).*max(rbktm(k,j))-a(j).*a(j)); % breakout area is cm^2 msld(j)=rhos*vst*area_sand(j); % in g/s else area_sand(j)=0; msld(j)=0; end msld_tot=sum(msld); end end for j=nst:1:nend for k=1:Tt if pwmid(k,j)<=pwc(k,j) mfld(j)=(pi*rhostar*max(vfmid(k,j)).*max(rbmx(k,j)).*max(rbktm(k,j)))*100*6.4516; % in g/s else mfld(j)=(rhof*pi*max(vfmid(k,j)).*a(j).*a(j))*100*6.4516; end end end % Fluid Density for j=nst:1:nend for k=1:Tt if pwmid(k,j)> pwc(k,j) rho(j)=rhof; %in g/cm3 elseif pwmid(k,j)<=pwc(k,j) && msld(j)mfld(j) rho(j)=rhof*((1-vol_sands(j)./vol_well(j))+(vol_sands(j)/vol_well(j))*rhos); %in g/cm3 end rhom=1000*sum(rho(j).*vol_well(j))/Volwel_ophl; % in kg/m3 end end model.param.set('rhom',rhom); % Fluid Viscosity vav=1/((nend-nst+1))*sum(max(vfmid(k,j))); %average velocity within the openhole section of the wellbore in m/s for j=nst:1:nend for k=1:Tt if pwmid(k,j)> pwc(k,j) muf(j)=mufw; elseif pwmid(k,j)<=pwc(k,j) && msld(j)2100 Cmu=0.09; elseif Nre<2100 Cmu=0; end model.param.set('Cmu',Cmu); %Rock Compressibility for j=nst:1:nend if index==1 cf(j)=97.32*10e-6/(1+79.8181*0.699993*por(j)); else cf(j)=0.8535/(1+2.202*10e6*1.075*por(j)); end end % Total Compressibility for j=nst:1:nend if M(j)==1 pbar(j)=0.5*(pint(j)+pwmidp(j));%average reservoir pressure cg(j)=1./pint(j); ct(j)=cf(j)+cg(j); else ct(j)=cf(j); end end hp=3.28084*[h1 h2 h3 h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15 h16 h17 h18 h19 h20 h21 h22 h23 h24 h25 h26 h27 h28];%true vertical thickness of the producing formation nstp=1; %number of the first producing layer nendp=10; %number of the last producing layer syms z for k=1:Tt for j=nstp:nendp %Gas Pseudo critical properties...according to Sutton Tpc(j)=169.2+349.5*gg(j)-74*gg(j).^2; Ppc(j)=756.8-131*gg(j)-3.6*gg(j).^2; %Pseudo Reduced properties Tr(j)=T(j)./Tpc(j); Pr(j)=pint(j)./Ppc(j); %Using Dranchuk and Abou-Kassem Z factor Correlation rhpr(j)=0.27*(Pr(j)./(Tr(j)*z)); %pseudo reduced density funz(j)=1+rhpr(j).*(0.3265-1.07./Tr(j)-0.5339./Tr(j).^3+0.01569./Tr(j).^4-0.05165./Tr(j).^5)... +(rhpr(j).^2).*(0.5475-0.7361./Tr(j)+0.1844./Tr(j).^2)-(rhpr(j).^5).*(-0.7361./Tr(j)+0.1844./Tr(j).^2)*0.1056... +0.6134*(1+0.721*rhpr(j).^2).*(rhpr(j).^2./Tr(j).^3).*exp(-0.721*rhpr(j).^2)-z; %Using Newton-Raphson Iterative Method in Estimating z-factor dfunz(j)=diff(funz(j),z); ddfunz(j)=diff(dfunz(j),z); % initializations funz0(j)=double(subs(funz(j),z,z0(j))); dfunz0(j)=double(subs(dfunz(j),z,z0(j))); ddfunz0(j)=double(subs(ddfunz(j),z,z0(j))); % Newton Raphson Convergence Criterion while(((funz0(j).*ddfunz0(j))/(dfunz0(j).^2))>1) z0=2*ones(1,nnp); funz0(j)=double(subs(funz(j),z,z0(j))); dfunz0(j)=double(subs(dfunz(j),z,z0(j))); ddfunz0(j)=double(subs(ddfunz(j),z,z0(j))); end z1(j)=z0(j)-(funz0(j)./dfunz0(j)); %check for accuracy while (abs(z1(j)-z0(j))>accz) z0(j)=z1(j); funz0(j)=double(subs(funz(j),z,z0(j))); dfunz0(j)=double(subs(dfunz(j),z,z0(j))); ddfunz0(j)=double(subs(ddfunz(j),z,z0(j))); z1(j)=z0(j)-(funz0(j)./dfunz0(j)); end bg(j)=(psc/tsc)*(z1(j).*T(j)./pint(j)); %average gas formation volume factor % Pseudo Pressure % Using Brown 1977 Correlation for Z- Factor n(j)=0.038-0.026*sqrt(Tr(j)./Tpc(j)); m(j)=0.51*(Tr(j)./Tpc(j)).^-4.133; %Lee, Gonzales, and Eakin 1966 for Gas Viscosity Ki(j)=((9.4+0.02*29*gg(j)).*Tr(j).^1.5)./(209+19*29*gg(j)+Tr(j)); Xi(j)=3.5+986./Tr(j)+0.29*gg(j); Yi(j)=2.4-0.2*Xi(j); %Pseudopressure Function psif=@(p) p./((1-m*(p/Pc)+n*(p/Pc).^2+0.0003*(p/Pc).^3).*(1e-4*Ki*exp(Xi*(29*g*p./((1-m*(p/Pc)+n*(p/Pc).^2+0.0003*(p/Pc).^3)*10.732*TR*62.4)).^Yi))); mpint(j)=2*integral(psif,pb,pint(j));%pseudo pressure of initial reservoir pressure mpwmidp(k,j)=2*integral(psif,pb,pwmidp(k,j));%pseudo pressure of wellbore flowing pressure of producing formations % Equivalent Radius for Inflow Due to Wellbore Enlargement Using % Equivalent Curved Surface Area Hypothesis hcp(j)=((rwp1(j)-rwp2(j))./(rwp1(j)+rwp2(j))).^2; rwpr(j)=0.5*0.0833333*(rwp1(j)+rwp2(j)).*(1+3*hcp(j)./(10+sqrt(4-3*hcp(j)))); % Using Srinivasa Ramanujan Approx. (in feet) %Applicable to Horizontal Wells td(j)=0.0002637*kh(j).*k*(24*Tt)./(por(j).*mur(j).*ct(j).*L(j)*0.5); %dimensionless time hd(j)=sqrt(kh(j)./kv(j)).*(2*hp(j)./L(j)); %dimensionless reservoir thickness zwd(j)=sqrt(kh(j)./kv(j)).*(2*zw(j)./L(j)); %dimensionless wellbore vertical position rwd(j)=rwpr(j)./L(j).*(1+sqrt(kh(j)./kv(j))); %dimensionless equivalent wellbore radius xmid(j)=0.5*(x1(j)+x2(j)); %Partial Penetration Skin Factor Pxyz(j)=(bL(j)./L(j)-1).*(log(hp(j)./rwpr(j))+0.25*log(kh(j)./kv(j))-log(180*zw(j)./hp(j))-1.84); if aw(j)/sqrt(kh(j))>=0.75*bL(j)/sqrt(kh(j))&& 0.75*bL(j)/sqrt(kh(j))>0.75*hp(j)/sqrt(kv(j)) %Wide Reservoir FL(j)=-(L(j)./2*bL(j)).*(0.145+log(L(j)./(2*bL(j)))-0.137*(L(j)./(2*bL(j))).^2); if (4*xmid(j)+L(j))/2*bL(j)<1 FLp(j)=-((4*xmid(j)+L(j))./2*bL(j)).*(0.145+log((4*xmid(j)+L(j))./2*bL(j))-0.137*((4*xmid(j)+L(j))./2*bL(j)).^2); else Flp(j)=(2-(4*xmid(j)+L(j))./2*bL(j)).*(0.145+log(2-(4*xmid(j)+L(j))./2*bL(j))-0.137*(2-(4*xmid(j)+L(j))./2*bL(j)).^2); end if (4*xmid(j)-L(j))/2*bL(j)<1 FLn(j)=-((4*xmid(j)-L(j))./2*bL(j)).*(0.145+log((4*xmid(j)-L(j))./2*bL(j))-0.137*((4*xmid(j)-L(j))./2*bL(j)).^2); else Fln(j)=(2-(4*xmid(j)-L(j))./2*bL(j)).*(0.145+log(2-(4*xmid(j)-L(j))./2*bL(j))-0.137*(2-(4*xmid(j)-L(j))./2*bL(j)).^2); end Pxy_prime(j)=2*bL(j).^2./(L(j).*hp(j)).*sqrt(kv(j)./kh(j)).*(FL(j)+0.5*(FLp(j)-FLn(j))); SR(j)=Pxyz(j)+Pxyz_prime(j); elseif bL(j)/sqrt(kh(j))>=1.33*aw(j)/sqrt(kh(j)) && 1.33*aw(j)/sqrt(kh(j))>hp(j)/sqrt(kv(j)) Py(j)=(6.28*bL(j).^2./(aw(j).*hp(j))).*sqrt(kv(j)./kh(j)).*((1/3-xmid(j)./bL(j)+(xmid(j)./bL(j)).^2)+L(j).*(L(j)./bL(j)-3)./(24*bL(j))); Pxy(j)=(bL(j)./L(j)-1).*((6.28*aw(j)./hp(j)).*sqrt(kv(j)./kh(j))).*(1/3-yw(j)./aw(j)+(yw(j)./aw(j)).^2); SR(j)=Pxyz(j)+Pxy(j)+Py(j); else SR(j)=0; end if inclp(j)=mpint(j) qsc(k,j)=0; end qinfl(k,j)=qsc(k,j).*bg(j)*3.27741e-4; %cubic meter per second else qinfl(k,j)=1.38009805e-6*(pint(j)-pwmidp(k,j))*kh(j)*hp(j)/(162.6*mur(j)*(log10(k*24*Tt)+log10(kh(j)/(por(j)*mur(j)*ct(j)*rwpr(j)*rwpr(j)))-3.23+ST(j))); %cubic meter per second if pwmidp(k,j)>=pint(j) qinfl(k,j)=0; end end else if td(j)<1 && td(j)25 % Early Time Radial Flow if M(j)==1 qsc(k,j)=(mpint(j)-mpwmidp(k,j)).*sqrt(kh(j).*kv(j)).*L(j)./(1638*T(j).*(log10(sqrt(kh(j).*kv(j)).*(k*24*Tt)./(por(j).*mur(j).*ct(j).*rwpr(j).^2))-3.2275+0.8686*SR(j)-2*log10(0.5*(nthroot(kh(j)./kv(j),4)+nthroot(kh(j)./kv(j),4))))); if mpwmidp(k,j)>=mpint(j) qsc(k,j)=0; end qinfl(k,j)=qsc(k,j).*bg(j)*3.27741e-4; else qinfl(k,j)=1.38009805e-6*(pint(j)-pwmidp(k,j)).*sqrt(kh(j).*kv(j)).*L(j)./(162.6*mur(j).*(log10(sqrt(kh(j).*kv(j)).*(k*24*Tt)./(por(j).*mur(j).*ct(j).*rwpr(j).^2))-3.2275+0.8686*SR(j)-2*log10(0.5*(nthroot(kh(j)./kv(j),4)+nthroot(kh(j)./kv(j),4))))); if pwmidp(k,j)>=pint(j) qinfl(k,j)=0; end end elseif (td(j)>zwd(j).^2 & td(j)=mpint(j) qsc(k,j)=0; end qinfl(k,j)=qsc(k,j).*bg(j)*3.27741e-4; else qinfl(k,j)=0.5*1.38009805e-6*(pint(j)-pwmidp(k,j)).*sqrt(kh(j).*kv(j)).*L(j)./(162.6*mur(j).*(log10(sqrt(kh(j).*kv(j)).*(k*24*Tt)./(por(j).*mur(j).*ct(j).*rwpr(j).^2))-3.2275+0.4243*SR(j)-log10((1+sqrt(kh(j)./kv(j))).*zw(j)./rwpr(j)))); if pwmidp(k,j)>=pint(j) qinfl(k,j)=0; end end % The rest of the flow periods are only applicable to no-flow % boundaries...top and bottom elseif td(j)>hd(j).^2 & td(j)>1 %Intermediate Linear Flow if G(j)==1 Sz(j) = -1.1513*sqrt(kh(j)./kv(j)).*(2*hp(j)./L(j)).*log10((pi*rwr(j)./hp(j)).*(1+sqrt(kv(j)./kh(j))).*sin(pi*zw(j)./hp(j))); %partial penetration pseudo skin in vertical direction if M(j)==1 qsc(k,j)=(mpint(j)-mpwmidp(k,j))./((81.9*z1(j).*T(j)./(L(j).*hp(j))).*sqrt(mur(j).*(k*24*Tt)./(por(j).*kh(j).*ct(j)))+1424*T(j)./(kh(j).*hp(j)).*Sz(j)+1424*T(j).*SR(j)./(L(j).*sqrt(kh(j).*kv(j)))); if mpwmidp(k,j)>=mpint(j) qsc(k,j)=0; end qinfl(k,j)=bg(j).*qsc(k,j)*3.27741e-4; else qinfl(k,j)=1.38009805e-6*(pint(j)-pwmidp(k,j))./((8.128./(L(j).*hp(j))).*sqrt(mur(j).*(k*24*Tt)./(por(j).*kh(j).*ct(j)))+141.2*mur(j)./(kh(j).*hp(j)).*Sz(j)+141.2*mur(j).*SR(j)./(L(j).*sqrt(kh(j).*kv(j)))); if pwmidp(k,j)>=pint(j) qinfl(k,j)=0; end end end else if (td(j)>1 & td(j)>hd(j).^2) if G(j)==1 Sz(j)=-1.1513*sqrt(kh(j)./kv(j)).*2*hp(j)./L(j).*log10((pi*rwpr(j)./hp(j)).*(1+sqrt(kv(j)./kh(j))).*sin(pi*zw(j)./hp(j)))-2*kh(j)./kv(j).*(hp(j)./L(j)).^2.*(1/3-zw(j)./hp(j)+(zw(j)./hp(j)).^2); if M(j)==1 qsc(k,j)=(mpint(j)-mpwmidp(k,j))./((1638*T(j)./(kh(j).*hp(j))).*(log10(kh(j).*(k*24*Tt)./(por(j).*mur(j).*ct(j).*(0.5*L(j)).^2))-2.5267)+1424*T(j).*Sz(j)./(kh(j).*hp(j))+1424*SR(j)./(L(j).*sqrt(kh(j).*kv(j)))); if mpwmidp(k,j)>=mpint(j) qsc(k,j)=0; end qinfl(k,j)=bg(j).*qsc(k,j)*3.27741e-4; else qinfl(k,j)=1.38009805e-6*(pint(j)-pwmidp(k,j))./((162.6*mur(j)./(kh(j).*hp(j))).*(log10(kh(j).*(k*24*Tt)./(por(j).*mur(j).*ct(j).*(0.5*L(j)).^2))-2.5267)+141.2*mur(j).*Sz(j)./(kh(j).*hp(j))+141.2*mur(j).*SR(j)./(L(j).*sqrt(kh(j).*kv(j)))); if pwmidp(k,j)>=pint(j) qinfl(k,j)=0; end end end end end end % Inflow Velocity...to be linked to COMSOL vinfl(k,j)=qinfl(k,j)./(2*0.0929*pi*rwpr(j)*L(j)); %in m/s end end %Linking the Flow Velocity with Comsol Deskstop model.param.set('vinfl1',vinfl(1,1)); model.param.set('vinfl2',vinfl(2,1)); model.param.set('vinfl3',vinfl(3,1)); model.param.set('vinfl4',vinfl(4,1)); model.param.set('vinfl5',vinfl(5,1)); model.param.set('vinfl6',vinfl(6,1)); model.param.set('vinfl7',vinfl(7,1)); model.param.set('vinfl8',vinfl(8,1)); model.param.set('vinfl9',vinfl(9,1)); model.param.set('vinfl10',vinfl(10,1)); model.param.set('vinfl11',vinfl(11,1)); model.param.set('vinfl12',vinfl(12,1)); model.param.set('vinfl13',vinfl(13,1)); model.param.set('vinfl14',vinfl(14,1)); model.param.set('vinfl15',vinfl(15,1)); model.param.set('vinfl16',vinfl(16,1)); model.param.set('vinfl17',vinfl(17,1)); model.param.set('vinfl18',vinfl(18,1)); %model.param.set('vinfl19',vinfl(19,1)); model.param.set('vinfl20',vinfl(20,1)); % model.param.set('vinfl21',vinfl(21,1)); model.param.set('vinfl22',vinfl(22,1)); model.param.set('vinfl23',vinfl(23,1)); model.param.set('vinfl24',vinfl(24,1)); model.param.set('vinfl25',vinfl(25,1)); % model.param.set('vinfl26',vinfl(26,1)); model.param.set('vinfl27',vinfl(27,1)); model.param.set('vinfl28',vinfl(28,1));