Estimation of in Situ Reserves Using Stochastic Finite Element
stochastic finite elements IN a complex system: Vibrating
NON STATIONARY fluid Flow, Mass TRANSPORT, heat Conduction.
Skender Osmani, Margarita .Qirko
Polytechnic University ,Tirana, ALBANIA
Abstract.
The paper presents a view on applications of stochastic finite elements(S.F.E.) R3 in one complex system, considering a differential equation with random parameters,
where,
x, y, z are the 3D spatial coordinates,. t- the processes time
is a process function( temperature, pressure etc)
Kij- the conductivity tensor.In general case it is:
( 12)
,d- capacity coefficient (function),g- mass coefficient (function),Q- density of the volume flux,
V –velocity vector.
This equation in specific conditions goes to:
- non stationary fluid flow in porous medium (or / and) ,
- mass transport (or /and)
- heat transfer (or/and)
- vibrating system etc.
The paper contains:
1-An approach of the mean value estimation of the parameter distributions, using S.F.E..
The next is defined as a block vi, with the random function Z(x), where xvi. is a random variable i.e. the value at point x determines the respective probability distribution p(x) .After the estimation of the mean value over the domain vi , is calculated by
zvi = (1)
2.-Development of the numerical model using S.F.E . applying a mixed algorithm at the , i.e. it is applied . the Galerkin,s approach not “as a whole” as it is often happened in the literature[3][15] , but partly , combining it with other numerical procedure as Runge-Kutta of the fourth order etc.
In this treatment, the initial and boundary conditions have been supposed to be treated specifically according by the given problem.
3.- Several simple examples as particular case of the mentioned equation.
4.Conclusions, the good things of the S.F.E in stationary flow, mass transport, heat transfer, vibrating , subsidence, waving, deformations, consolidations, earthquakes and other phenomenon’s.
1.Introduction.
Some processes are complex, as non stationary flow, mass transport, heat transfer, vibrating , subsidence, waving, deformations, consolidations, earthquakes and other geodynamics , reservoir engineering [1][14][16]
Their complexity depend on their random mechanical, physical, chemical etc parameters, relations between them, their risk state etc. .Several processes could be deterministic, other ones stochastic, several could be studded “separately and independently” as the simple heat conduction, one phase fluid flow in porous medium etc.
But in some circumstances the situation seems to be more complex as it the case of tsunami which is a large destructive sea wave generated by geodynamics phenomenon’s as earthquake, subsidence , volcanic eruption..
Of course in really conditions the above phenomenon’s are often characterized by the their risk and uncertainty of initially and boundary conditions etc, so the situation could be more aggravated.
Taking into account the above considerations as well as to simplify as much as possible the treatment of , we choose a differential equation as “ a general case “ for SFE applying.
2.Stochastic Finite Element.
2.1Mean value estimation
Let’s consider: a zone V R3 and a random function Z(x) , xV; the zone V is partitioned into blocks vi by a parallelepiped grid :
V= vI (2)
where vi is a parallelepiped element with eight nodes. In each node Z(x) is known i.e the probability density on its point support as it is shown in Fig(1) .
Fig. 1. Parallelpiped element.
Let, suppose ,it is required 8:
-the distribution p.d.f in whatever point xV
-the estimation of the mean value
zvi = 1/vi Z(x)dx (3)
over the domain vi .
The stochastic finite element is defined as a block vi, with the random function Z(x), where xvi. is a random variable i.e. the value at point x determines the respective probability distribution p(x).
Let us consider a reference element wi in the coordinate system s1 s2 s 3 2. If we choose an incomplete base:
P(s)=1 s1 s2 s3 s1s2 s2s3 s3s1 s1s2s3 (4)
then the function Z(x) could be presented as a linear combination[ ] :
Z(x)= Z( s1 s2 s 3 )= < P(s) > [ P8 ] –1 { Zs 8 } = < N (s) >{ Zs 8 } (5)
where [ P8 ] –1 is the matrix , whose elements are the polynomials base values at the nodes,{ Zs 8 } is the vector of the distributions of the nodes, while N (s) is the vector of the shape functions
Ni, i= 1,8
To calculate the mean value zvi = 1/vi Z(x)dx , we consider the deterministic transformation :
Xi (s)= < Ni (s) > < xi8 i=1,3 ( 6)
Therefore
Zvi =Z (x1(s1 s2 s3 ), x2(s1 s2 s3 ), x3(s1 s2 s3 )) det J ds1 ds2 ds3=
< N (s) >{ Zs 8 } det J ds1 ds2 ds3 = HiZi(x) (7)
where J is the Jacobian of the transformation , Hi are the distribution weights depending only on the node coordinates.and . In other words they make the weighted average of the given distributions at the nodes. Furthermore , if we consider the expectation:
E{Zv(x)}=E((HiZi(x)dx (8)
under the hypothesis:
E (Zi(x) )= mi m = const , i=1,8 (9)
then
E{ Zv(x)}= Hi E (Zi(x) ) = Him=m , i=1,8 (10)
Thus,the stochastic estimator of chose finite element is a linear interpolator, related to the distributions given at its nodes.
A small example is presented to illustrate the procedure of mean value estimation[6]. The distributions of temperature placed at the nodes of a quadrilateral are shown in figure (2)(3). Applying the code basing on the above algorithm, the mean distribution and cumulative function is shown in figure (4)(5)
Fig 2 Node density distribution Fig 3 Node cumulative function
Fig 4 Mean value density distribution Fig5 Mean value cumulative function
The temperature mean value of around 140C, has the highest probability, i.e the most probable energy volume.
2.2Numerical model
Let’s briefly consider some differential equations[12] in a zone c R3
1-Heat equation
Q (11)
where,
x, y, z are the 3D spatial coordinates.
T-the temperature function: T=T(x, y, y, z, t)
Kij- the conductivity tensor.In general case it is:
( 12)
t- the processes time,d- capacity coefficient (function),g- mass coefficient (function),Q- density of the volume flux
2-Solute transport equation
(13)
(13)
c- concentration in the liquid phase as mass of contaminant per unit volume of solution
Dx, Dy, Dz, - directional hydrodynamic dispersion coefficient.,R- retardation factor depending on other factors [15][16] ,k- the overall first order decay rate,Vx, Vy, Vz, -directional (seepage) velocity components
3.Vibrating (system) equation:
(14)
m-mass of the system,d- damping coefficient (function),k- elasticity coefficient (function),f- external force[14]
As the applied SFE mechanism is analogues to above equations, below we are considering another equation in a given zone c R3
(15)
The proposed equation is a mixed one, composed by parabolic, hyperbolic etc equations [17], so in specific conditions it could behave and work like one of them. In this context this equation could be considered as “a general” one.
We assume the equation has its parameters like random variables in whatever point M.We divide the zone c R3 into parallelepiped elements, one of which is shown in figure (1) We shall now solve the differential equation in one SFE, afterward the same procedure will be repeated for the whole structure of SFE.
In the SFE nodes are known the respective distribution, so we can suppose the mean value of those distributions has been obtained, even in which node we know the part of it. This remains to find only the solution already as if the parameters of equation are deterministic values.
In this context, we can first apply the Galerkin’s approach, afterwards the Runge- Kutta of the fourth order.For the parallelepiped element we have:
(16)
where <N> the vector of the shape functions, while {}the vector of the model values of .According to Galerkin’s approach, we multiply the two sides of the equation by Nje, j=1, 2, -----8. Thus,
(17)
where (nabla operator) (18)
Applying the Green theorem to the first term of the equation
(19)
and carrying out the algebraic operations, we obtain:
, i, j = 1,8 (20)
If we note:
(21)
then
(22)
Let’s assume the directional velocity components are as following [15]:
(23)
where:
h- function (potentrometric head) which is the solution of the differential equation:.
(24)
where Cx, Cy, Cz, KLxx, KLzz, are the parameters, in several cases random variables.[8][9]. Again, applying the Galerkin’s method and Green theorem, the equation becomes:
(25)
where i,j =1.8 (26)
(27)
{Fh}- the vector resulting from the Green theorem. In the other hand,
(28)
According to the Galerkin’s approach we have:
(29)
The above equations could be written in matrix form as:
(30)
where,, are the coefficient matrices, while the vectors ,,are column vectors containing boundary conditions.
Taking into account all system of equations as well as the boundary conditions we can find the vectors ,,. This, the velocities Vx, Vy, Vz be considered as the known in the differential equation’s
So in these circumstances the velocities, Vx, Vy, Vz will be considered known. Being so, further we will again simplify the idea, considering only the differential equations system.
(31)
Supposing that system has been written for every element and assembling all these systems into one for the whole structure we have:
(32)
Already the vector {} depends only on the variable t:
(33)
where (34)
Noting: (35)
then
(36)
These considerations lead to the problem of finding the solution of a system, which contains only the differential equations of the first order. The total number of the system equations is equal to the number of the finite element nodes in the whole structure[12].On this circumstances knowing the specific initial and boundry conditions we can find the respective solution.
Basing on the mentioned algorithm with three dimension SFE it is compiled a C++ code for the general equation(15).
Below, it’s applications are presented some illustrative examples in which the parameters of the equations are considered as deterministic ones. That is , because it is assumed the respective given distributions are already estimated according the above mean value approach. Thus , we will simply work with deterministic parameters, as well as with deterministic initial and boundary conditions.
3.Illustrative examples.
3.1.Vibrating system:
In one SFE structure are considered only two nodes of a vibrating system[13][14] , characterized of the following differential equations:
d2u/dx2-3du/dx +4u=6sin(2t)
d2v/dx2-3dv/dx +2v=exp(2t)*sin(t) (37)
The initial conditions of the functions and their derivatives are :
u0=1 , u0’=0, (38)
v0=-0.4, v0’=-0.6
The solution is asked at the process time t[0, 1].In order to obtain it,we take into account the Runge Kutta approach and transform the system as bellow:
du1/dt = u2 = f1(t,u1,u2,u3.u4) (39)
du2/dt = 6sin(2t)+3u2-4u1 = f2 (t,u1,u2,u3.u4)
du3/dt = u4 = f3(t,u1,u2,u3.u4)
du2/dt = exp(2t)*sin(t)+2u4-2u3= f4 (t,u1,u2,u3.u4)
For different time :
t=0.1, 0.2, 0.3, 0.4, 0.5
(dividing the given segment into elementary ones with n=4000),are obtained the following results:
U[0.1]=0.9801 U[0.2]= 0.9212 U[0.3]= 0.8256 U[0.4]= 0.6971 U[0.5]= 0.5408
V[0.1]=-0.4617 V[0.2]=-0.5256 V[0.3]=-0.5886 V[0.4]=-0.6466 V[0.5]=-0.6936
Taking into account the exact solution of the given system:
u=cos(2*t);
v=0.2(exp(2t)*(sin(t)-2*cos(t))
the calculation error of the solution results generally smaller than 0.0001.
It is shown, the above algorithm could be used for the resolution of differential equations systems
independently if they are providing from finite element structures or others ones. In other words
this algorithm “does not belong” only to SFE but also to others “ mechanisms”.
3.2.Non stationary heat transfer:
As a application of heat transfer let’s consider the plane problem :
cp*(d 2U/dx2+k*d2U/dy2)+cd*dU/dt +cg*d2U/dt2 = Q(x,y,t) (39 )
cp=1,cd=1,cg=1;
in one SFE face in which three nodes having their coordinates as:
1 x=0,y=1; 2 x=0,y=0; 3. x=1,y=0 (39.1)
Initial conditions :
U1-Initial function U2- Initial derivative
U1[1]=1 U2[1]=-3 (39.2)
U1[2]=0 U2[2]=-3
U1[3]=1 U2[3]=-3
Boundary conditions: 1 x=0,y=1 ,U=1, du/dt=-3 ; t>0;
3 x=1,y=0 ,U=1, du/dt=-3 ; t>0; (39.3)
The question is to calculate the function U in the node 2 for the time t: 0<t<=10 .Applying the mention C++ code there are obtained:
Vector U1 Vector U2
t=0
U1[1]=1 U2[1]=-3
U1[2]=0 U2[2]=-3
U1[3]=1 U2[3]=-3
t=2
U1[1]=-5.0000003 U2[1]=-3
U1[2]=-4.4563241 U2[2]=-3.6538434 (39.4)
U1[3]=-5.0000003 U2[3]=-3
t=5
U1[1]=-14.000001 U2[1]=-3
U1[2]=-13.928011 U2[2]=-3.3199201
U1[3]=-14.000001 U2[3]=-3
t=10
U1[1]=-29.000001 U2[1]=-3
U1[2]=-28.866354 U2[2]=-3.0382929
U1[3]=-29.000001 U2[3]=-3
The exact solution of this problem is:
U=x2+y2-3*t.
Comparing theoretical results to those obtained , the error of heat transfer gets to be smaller for t>10.
Further , lets, suppose at the two points 1, 3 are located two sources ( say exploration, exploitation wells, thermal, renewable energy sources etc ) while we are interested to design a project for another well. The point is: basing on the information of two sources Nr 1,3 of our example ,where are going to be the best (optimal) [4][10] position of the source Nr2. from economical view ?.
SFE model is an easy and simple alternative solution, without asking for the supplementary information in the given zone related to the exploration, exploitation, well testing[1] etc.It don’t need for heavy finite element structure , starting form natural boundary condition etc Also it works step by step , from known nodes to another unknown basing on the source equilibrium,. Even more the SFE model could and should be complementary to other ones as geostatistical , risk analysis etc.
3.3 Heat general equation..
Now , in one reference [ ] SFE , let’s consider the general equation of the heat transfer :
cp1*d2u/dx2+cp2*d2u/dy2+cp3*d2u/dz2 + cg*d2u/dt2+cd*du/dt = Q(x,y,z,t)
: cp1= cp2= cp3=cd=cg=1;
.
Element node coordinates are as following :
M1( -1,-1,-1), M2(1,-1,-1), M3( 1,1,-1),M4(-1,1,-1), M5( -1,-1,1), M6(1,-1,1), M7( 1,1,1), M8(-1,1,1);
-Initial conditions :
In the code , initial conditions in all nodes 1,2,3,4,5,6,7,8 are generated from the function :
U=(x2+y2+z2)exp(-t), t=0.
-Boundary conditions
a-Also, the function values in the boundary nodes 1,2,3,4 are generated from the function :
U=(x2+y2+z2)exp(-t) ; .t>0.
The question is to find the function in the nodes 5,6,7,8.
The respective results (given by the mentioned C++ code) are as following:
Vektori U1 Vektori U2
t=0.1
U1[1]=2.7145123 U2[1]=-2.7145123
U1[2]=2.7145123 U2[2]=-2.7145123
U1[3]=2.7145123 U2[3]=-2.7145123
U1[4]=2.7145123 U2[4]=-2.7145123
U1[5]=2.7372444 U2[5]=-2.1152337
U1[6]=2.7355745 U2[6]=-2.1417675
U1[7]=2.7355828 U2[7]=-2.1414897
U1[8]=2.7348073 U2[8]=-2.1599109
t=0.30
U1[1]=2.2224547 U2[1]=-2.2224547
U1[2]=2.2224547 U2[2]=-2.2224547
U1[3]=2.2224547 U2[3]=-2.2224547
U1[4]=2.2224547 U2[4]=-2.2224547
U1[5]=2.4427445 U2[5]=-0.8000741
U1[6]=2.4322803 U2[6]=-0.84876978
U1[7]=2.43223 U2[7]=-0.85023367
U1[8]=2.4262745 U2[8]=-0.88062537
The exact solution of the above differential equation is U=(x2+y2+z2)exp(-t), so the solution has a spheric symmetry .This fact is accomplished only for t=0.1 , while for t=0.3 that is not.
However the result are acceptable, in the case analyzed it is possible to increase the accuracy of the solutions intervening at the geometrical, time dimensions etc.
As in the second example , the solution of two sources unknown (sides 5-6,7-8) is provided only on base of two other known sources (sides 1-2,3-4) and not more.
4.Conclusions.
-The study of complex phenomenon’s could be done using SFE. Their applications are essential in the geodynamics, energy[8] and water resources, reservoir engineering( exploration , exploitation, oil, gas water flow[6], well testing[1] ) , hydrology, environment, geostatistics, risk analysis[11] etc
- SFE algorithm could be used not only in the solution of the mentioned equations ,but also in geostatistics (Kriking estimation),plane and space stress, bending etc [12],in the local estimation of different proprieties, solving the inverse local problems etc .The geostatistics methods could be simultaneously used with other dynamic ones for parameter estimation , reserve calculation[7] [9] etc .
- SFE could be coded using Object – Oriented Languages as C++[5], independently of (excessive) data, complex and difficulties programs, which in the future will need to be modified , extended etc as it was the case analyzed etc.
- In the view of resources management [10], the SFE approach could and should be used as an important strategy in the exploration, exploitation, using , monitoring, reserve estimation, optimal placement and location etc of energy, oil, water etc. resources in different engineering projects, shortly say to determine and decide of optimally economical solution, the best of the contingency alternatives .
References
[1] Dake L.P. Fundamentals of Reservoir Engineering Nineteenth edition, Elsevier 2002,UK
[ 2]Dhat.G , Touzot G. The finite element method displyed J.Willey & Sons, Paris1984.,France
[ 3] Future Groundwater Resources at Risk. Proceedings of the third International Conference,
Lisbon 2001, Portugal.
[ 4]Gregoire A. Analyse numerique et optimization, Ec. Polytechnique .Paris 2003, France
[ 5 ]Gaddis T. Krupnow B.Starting out with C++.( fourth edition) Scott /Jones 2004,USA
[ 6]Osmani.S.Doracaj.M. Gottardi G. A three phase, threedimension finite element.
Bolletinno Asocc.Mineraria Subalpina Nr.1,2 Torino 1990, Italy.
7 Osmani.S. Finite Element on Geostatistique Calculations. In A.G. Fabri et Royer , 3 rd International Codata Conference in Geomathematics and Geostatistics Sci. de la terre , Ser Inf. Nancy ,1994:32 France
8 Osmani.S. Energy distribution estimation using Stochastic Finite Element
Renewable Energy . Elsewier 25 (2002) , London, UK.
[ 9] Osmani S. Qirko. M. Stochastic Finite Element for petroleum reservoir properties estimation
3rd International Conference of Applied Geophysics, Cairo, 2006. Egypt
10 Osmani.S..Osmani F. Consideration on the Strategy , Policy and Management of Resources.
Paper accepted to be presented at the International Conference : Survival and Sustainability
19-24 February ,2007, Nicosia., Northern Cyprus.
[11]Osmani .S,Aleti R.Geostatistics,Finite element and Risk Analysis in Environmental Resources Estimation :International Conference on..Environmental Problems of Mediterranean Region.,Nicosia 2002., Northern Cyprus.
[ 12 ]Osmani S. Finite Element (in Albanian) Vol 1,2 Sh.B.L.U , Tirana 1988., Albania
[ 13] Qirko M.Numerical Analysis(in Albanian) Sh.B.L.U , Tirana 2004., Albania
[ 14 ]Verruijt A.Soil Dynamics . Delft University 2005.The Netherland.
[ 15 ] Sherif M.M Singh.V.P,A.Rashed. Proceedings of the International Conference on Water Resources
Vol 3,4,2002. Kuwait
[16 ] Fundmantals of transport phenomena in porous media IAHR.Elsevier 1972 London
[17 ] Colton D. Partial differential equations , New York. 1988,USA
.
1