Download Computer Project 4: Approximating the Integral of Ei(x) for Geosystems Engineering and more Study Guides, Projects, Research Engineering in PDF only on Docsity!
SOLUTION SET
Computer Project No. 4
February 12, 2004 Due on February 17, 2004, 9:30 AM
PGE310 (Unique No. 17235) Spring Semester 2004 Formulation and Solution of Geosystems Engineering Problems
Instructor: Carlos T. Verdín
MORE ON INPUT/OUTPUT COMMANDS AND CONSTRUCTION
OF FUNCTIONS WITH MATLAB
DESCRIPTION:
This computer project is intended to expand the programming exercises with for loops and flow control operations used for the solution of Computer Project No.
- In addition, the project introduces the notion and use of function subprograms as a means to improve your programming efficiency with Matlab.
YOUR REPORT SHOULD BE CLEAN, NEAT, AND WELL ORGANIZED; IT SHOULD INCLUDE CLEAR AND COMPLETE TECHNICAL DESCRIPTIONS WHEREVER NECESSARY. ALL FIGURES AND TABLES SHOULD BE LABELED AND PROPERLY ANNOTATED WITH A CAPTION. LOOSE FIGURES ARE NOT SELF-EXPLANATORY. POINTS WILL BE DEDUCTED FROM PROJECTS THAT DOES NOT ADHERE TO THESE PRESENTATION RULES.
Problem No. 1
Hydrocarbon producing wells are often shut in for a number of reasons, some of them having to do with regular maintenance operations such as workovers. On occasion, wells are also shut-in in order to generate reservoir pressure transients. The latter are used to infer formation properties such as permeability, porosity, spatial connectivity, fractures, etc. There are two cycles for reservoir pressure transients, referred to as drawdown and buildup. When the reservoir is being produced, pore pressure decreases because of the reduction in fluid support. Conversely, when the reservoir is shut-in reservoir pressure increases because of the equilibration of pressure support by pore fluids. Drawdown
pressure transients occur when the reservoir is shut in, while buildup pressure transients occur when the reservoir is put back in production mode. Both cycles of pressure transients are used to infer reservoir petrophysical properties. Pressure transients are measured either in the same well or else in a different well within the same reservoir.
When a well is producing reservoir fluids at a constant rate of qB and is suddenly
shut-in, the induced pressure drawdown is given by
( )
948 2 e ( , ) 70.6 where for 0
x^ u t i
qB c r p r t p Ei Ei x du x kh kt u
μ φμ −∞
∫ ,^ (1)
[One can show that Ei (^) ( − x ) (^) = − Ei x ( ).]
where (the values assigned to the variables below are examples): q is the production rate, q = 600 STB day / B is the oil formation volume factor, B =1.2 RB STB / p is the pressure in psi, pi = 5000 psi r is the radial location measured from the borehole axis, r = 10 ft t is time measured in hour, 0.01 < t < 1000 hr φ is porosity measured in fractions, φ =20 % k is permeability measured in md, k = 200 md
μ is the viscosity in cp, μ =1.5 cp
h is the formation thickness, h = 80 ft c t is the total compressibility, 30 106 1 ct psi = × −^ −
Equation (1) above assumes an infinite well and a homogeneous and isotropic reservoir. It also assumes that there is only one fluid phase involved in the production.
The function Ei x ( (^) ) in equation (1) does not have an analytical solution. It has to
be approximated with elementary functions and/or polynomial expansions. One possible approximation is obtained by means of the Taylor series expansion, namely,
( ) ( ) 1 (^ )
ln for 0 where =0.5772156649 is Euler's constant !
k
k
x Ei x x x k k
+∞
=
= + + (^) ∑ >
In practice, the above function has to be truncated to an order n^ , i.e.,
( ) ( ) 1 (^ )
ln !
n^ k
k
x Ei x x k k
=
≈ + + (^) ∑
sum =sum +term; end end y=sumsign(x); end*
The program is written in a different Matlab file than the function:
disp('') n = input('Enter the order of the truncation, n = ');
while n<0 | floor(n)~=n disp('n is a positive integer') disp(' ') n=input('Please re-enter n = ') end
% see function expint310.m
x = logspace(-1,0,101);
for i = 1:length(x) Ei(i) = expint310(x(i),n); end
figure(1); plot(x,Ei,'b'); title('Plot of the approximation of Ei(x) vs x'); xlabel('x'); ylabel('Ei(x)');
(3) Verify the accuracy and reliability of your Matlab function for Ei x ( ) by
comparing it to Matlab’s internal functions. Hint : to make use of Matlab’s utilities, use the identity Ei x ( (^) )= real(-expint(-x)) , where expint(x) is the internal Matlan function. The latter expression is ONLY valid for positive values of x. Compare the results of your Matlab function against those of the Matlab’s internal function as a function of x for x between 0.1 and 1. Quantify the errors of your Matlab function, if any.
% Comparing built-in Ei(x) and approximation
for i = 1:length(x) Ei_matlab(i) = real(-expint(-x(i))); end
relative_error = abs((Ei-Ei_matlab)./Ei_matlab)100;*
figure(2); plot(x,relative_error,'b'); title('Plot of the relative error between Ei(x) and its approximation vs x'); xlabel('x'); ylabel('Relative error, %');
Obviously, the error gets closer to zero when the order of approximation increases. The maximum relative errors of the approximation of Ei x ( (^) ) for n = 2 , n = 5 , n = 10 , n = 20 and n = 50 are respectively 80%, 0.017%, 1.3x10- (^7) %, 4.2x10 -8 (^) %, 4.2x10 -8 (^). The approximation seems to be acceptable for n > 5_. A peak is seen at x_ ≈ 0.3715 no matter what order of approximation was taken. The simulation for n = 5 and n = 10 show relatively high error for x → 1_. This is due to the fact that the approximation obtained from the Taylor series expansion is taken around x_ 0 (^) = 0_. Thus its approximation is best when x_ → 0_._
(4) Making use of your Matlab function for Ei x ( (^) ), compute and plot p ( , ) r t vs
log (^) ( ) t for three value of permeability, namely, k =50 md, k =400 md, and k =1000 md. Plot the three results as a function of time on the same graph with different colors. Make use of the time vector constructed in item (1) to construct your plots. Pay attention to units of pressure and time and use properly label the axes of your plot with the corresponding units. What can you conclude about the behavior of pressure as a function of permeability?
% Comparing built-in Ei(x) and approximation
q = 600; B = 1.2;
pi = 5000; r = 10; phi = 0.2; mu = 1.5; h = 80; ct = 0.00003; for i = 1:length(t) p_50(i)=pi+70.6qBmu/(50h)expint310(- 948phimuctr^2/(50t(i)),n); p_400(i)=pi+70.6qBmu/(400h)expint310(- 948phimuctr^2/(400t(i)),n); p_1000(i)=pi+70.6qBmu/(1000h)expint310(- 948phimuctr^2/(1000t(i)),n); end figure(3); plot(t,p_50,'b',t,p_400,'k--',t,p_1000,'m-.') title('Build-up plot of p vs t at different permeabilities for n=50'); xlabel('x'); ylabel('Relative error, %'); legend('k = 50md','k = 400md','k = 1000md',4);**
The lower k , the faster the pressure increases.
