A problem with Stehfest numerical laplace inverse algorithm

1 vue (au cours des 30 derniers jours)
edu
edu le 17 Mar 2014
Réponse apportée : edu le 17 Mar 2014
Hi. I've downloaded Stehfest numerical laplace inverse algorithm from this website and applied it to my own problem (i.e. Well testing in petroleum). The procedure is to use Stehfest algorithm in order to bring back pressure derivative response of the well from Laplace media to time domain. To do this, some additive parameters are included in the Laplace solution which include: Skin, dimensionless wellbore storage and etc. The Stehfest algorithm applies successfully for skin values greater than zero but when I use negative skin values a problem emerges. The problem is that using negative skin values in the laplace inverse equation , resulted pressure derivative response in time domain would be negative! The negative pressure derivative response has no physical meaning. I've checked Stehfest algorithm by some simple function which have exact analytical solution in time domain and the algorithm works well for those functions. Here is the form of laplace equation which should be converted into time domain. Your help is really appreciated.
function f=fpss(p)
% Finite Acting NFR
% PSS IPF
global a1 a2 a3 a4 a5
o=a1; lam=a2; Cd=a3; S=a4; red=a5;
% F(s) function=g
g=(o*(1-o)*p+lam)/((1-o)*p+lam);
a=sqrt(p*g);
b=a*red;
X=besseli(1,a)*besselk(1,b)-besseli(1,b)*besselk(1,a);
Y=besseli(0,a)*besselk(1,b)+besseli(1,b)*besselk(0,a);
%Pressure solution in laplace media
f1=(besseli(0,a)*besselk(1,b)+besseli(1,b)*besselk(0,a)...
+(S*a*(besseli(1,b)*besselk(1,a)-besseli(1,a)*besselk(1,b))))...
/((Cd*p^(2)*(Y-S*a*X))-p*a*X);
% Pressure derivative solution in laplace media
f=f1*p;

Réponses (1)

edu
edu le 17 Mar 2014
Actually I think the problem maybe related with double precision used by Matlab to do calculations and resulted round-off error. It may be proposed to use VPA , but this increases CPU time dramatically.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by