/[escript]/trunk/finley/test/python/coalgas.py
ViewVC logotype

Annotation of /trunk/finley/test/python/coalgas.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3484 - (hide annotations)
Thu Mar 24 05:42:07 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 19464 byte(s)
some work on coal gas
1 gross 3458 #######################################################
2     #
3     # Copyright (c) 2003-2010 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6     #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10     #
11     ########################################################
12     """
13     Gas in Coal Seam (fully coupled version)
14     """
15     __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16     Earth Systems Science Computational Center (ESSCC)
17     http://www.uq.edu.au/esscc
18     Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21     __url__="https://launchpad.net/escript-finley"
22    
23 gross 3484 from esys.escript.linearPDEs import LinearPDE
24 gross 3458 from esys.escript import unitsSI as U
25 gross 3484 from esys.escript import sqrt, log
26     import math
27 gross 3458
28 gross 3478 class MaterialProperty(object):
29     """
30     generic class for material properties depending on one (or several) status
31     variable(s)
32     """
33     def __init__(self, *args, **kwargs):
34     """
35     set up
36     """
37     pass
38     def __call__(self, *args, **kwargs):
39     """
40     return value of material property for given status
41     """
42     return self.getValue( *args, **kwargs)
43    
44    
45     def getValue(self, *args, **kwargs):
46     """
47     return value of material property for given status
48    
49     :remark: needs to be overwritten
50     """
51     raise NotImplementedError
52 gross 3458
53 gross 3478 class MaterialPropertyWithDifferential(MaterialProperty):
54     """
55     generic class for material properties depending on one (or several) status
56     variable(s) where the derivative of the property with respect to the
57     status variables is available
58     """
59     def getValueDifferential(self, *args, **kwargs):
60     """
61     returns the value of the derivative of material property for given status variable
62    
63     :remark: needs to be overwritten
64     """
65     raise NotImplementedError
66 gross 3484
67     class Porosity(MaterialPropertyWithDifferential):
68     """
69     defines porosity phi as function of pressure
70    
71     phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
72    
73     phi_p_ref is claculted from the initial porosity phi_0 at pressure p_0
74     """
75     def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
76     """
77     """
78     self.phi_p_ref=1 # will be overwritten later
79     self.p_ref=p_ref
80     self.C=C
81     # reset phi_p_ref to get phi(p_0)=phi_0
82     self.phi_p_ref=phi_0/self.getValue(p_0)
83 gross 3478
84 gross 3484 def getValue(self, p):
85     """
86     returns the porosity for given pressure p
87     """
88     X= self.C * ( p - self.p_ref )
89     return self.phi_p_ref * (1. + X * (1. + X/2 ) )
90    
91     def getValueDifferential(self, p):
92     """
93     returns the porosity for given pressure p
94     """
95     X= self.C * ( p - self.p_ref )
96     return self.phi_p_ref * self.C * (1. + X )
97    
98    
99     class WaterDensity(MaterialPropertyWithDifferential):
100     """
101     set water density as
102    
103     rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
104    
105     with rho_ref =rho_s/B_ref * gravity
106     """
107     def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1., C = 0./U.bar, rho_s= 998.2 * U.kg/U.m**3):
108     self.rho_ref = rho_s/B_ref * gravity
109     self.C=C
110     self.p_ref=p_ref
111    
112     def getValue(self, p):
113     """
114     returns the density for given pressure p
115     """
116     X= self.C * ( p - self.p_ref )
117     return self.rho_ref * (1+ X * (1+ X/2) )
118    
119     def getValueDifferential(self, p):
120     """
121     """
122     X= self.C * ( p - self.p_ref )
123     return self.rho_ref * self.C * (1+ X)
124    
125     class WaterViscosity(MaterialProperty):
126     """
127     defines viscosity of water
128    
129     mu=mu_ref /(1 + X + X*X/2)
130    
131     with X=-C*(p-p_ref)
132     """
133    
134     def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
135     """
136     set up
137     """
138     self.mu_ref = mu_ref
139     self.p_ref = p_ref
140     self.C=C
141     def getValue(self, p):
142     """
143     returns the viscosity for given pressure p
144     """
145     X= -self.C * ( p - self.p_ref )
146     return self.mu_ref/(1+ X * (1+ X/2) )
147    
148     class GasDensity(MaterialPropertyWithDifferential):
149     """
150     set water density as
151    
152     rho = gravity * rho_air_s /B(p)
153    
154     where B is given by an interpolation table
155     """
156     def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
157     self.rho_ref =rho_air_s * gravity
158     self.tab=InterpolationTable(x=p, y=B)
159    
160     def getValue(self, p):
161     """
162     returns the density for given pressure p
163     """
164     return self.rho_ref/self.tab.getValue(p)
165    
166     def getValueDifferential(self, p):
167     """
168     """
169     B = self.tab.getValue(p)
170     dBdp = self.tab.getValueDifferential(p)
171     return -self.rho_ref * dBdp/(B * B)
172    
173     class InterpolationTable(MaterialPropertyWithDifferential):
174 gross 3478 """
175     a simple 1D interpolation table for escript Data object with non-equidistant nodes
176     """
177     def __init__(self,x=[], y=[], obeyBounds=True):
178     """
179     set up interpolation table. obeyBounds is set an exception is thrown if
180     the interpolation argument is below min(x) or larger than max(x). Otherwise
181     the value for x is set to y[0] for
182     """
183 gross 3484 MaterialPropertyWithDifferential.__init__(self)
184 gross 3478 if len(x) != len(y):
185     raise ValueError,"length of interpolation nodes and value lists need to be identical."
186 gross 3484 if len(x) < 1 :
187 gross 3478 raise ValueError,"length of interpolation nodes a list needs to at least one."
188    
189     x_ref=x[0]
190     for i in range(1,len(x)):
191     if x_ref >= x[i]:
192     raise ValueError,"interpolation nodes need to be given in increasing order."
193     x_ref=x[i]
194     self.__x = x
195     self.__y = y
196     self.__obeyBounds = obeyBounds
197 gross 3484
198 gross 3478 def getValue(self, x):
199     """
200     returns the interpolated values of x
201     """
202     x0=self.__x[0]
203     m0=whereNegative(x-x0)
204     if self.__obeyBounds:
205     if sup(m0) > 0:
206     raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
207     out=self.__x[0]
208     for i in range(1,len(self.__x)):
209     x1=self.__x[i]
210     z=(y[i]-y[i-1])/(x[i]-x[i-1]) * (x-x[i-1]) + y[i-1]
211     out = out * m0 + z * (1-m0)
212     m0=whereNegative(x-x[i])
213    
214     if self.__obeyBounds:
215     if inf(m0) < 1:
216     raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
217     else:
218     out = out * m0 + y[-1] * (1-m0)
219     return out
220 gross 3458
221 gross 3484 def getValueDifferential(self, x):
222     x0=self.__x[0]
223     m0=whereNegative(x-x0)
224     if self.__obeyBounds:
225     if sup(m0) > 0:
226     raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
227     out=0.
228     for i in range(1,len(self.__x)):
229     x1=self.__x[i]
230     z=(y[i]-y[i-1])/(x[i]-x[i-1])
231     out = out * m0 + z * (1-m0)
232     m0=whereNegative(x-x[i])
233    
234     if self.__obeyBounds:
235     if inf(m0) < 1:
236     raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
237     else:
238     out = out * m0
239     return out
240    
241    
242 gross 3478 class Well(object):
243     """
244     generic well
245    
246     :var WATER: phase identifier for water well
247     :var GAS: phase identifier for gas well
248     """
249     WATER="Water"
250     GAS="Gas"
251 gross 3484 def __init__(self, Q=0., schedule=[0.], phase="Water", *args, **kwargs):
252 gross 3478 """
253     set up well
254     """
255 gross 3484 self.__schedule=schedule
256     self.__phase=phase
257     self.__Q = Q
258    
259 gross 3478 def getGeometry(self):
260     """
261     returns a function representing the geometry of the well.
262     typically a Gaussian profile around the location of the well.
263    
264     :note: needs to be overwritten
265     """
266     raise NotImplementedError
267 gross 3458
268 gross 3478 def getProductivityIndex(self):
269     """
270     returns the productivity index of the well.
271     typically a Gaussian profile around the location of the well.
272    
273     :note: needs to be overwritten
274     """
275     raise NotImplementedError
276    
277 gross 3484 def getFlowRate(self, t):
278 gross 3478 """
279     returns the flow rate
280     """
281 gross 3484 if self.isOpen(t):
282     return self.__Q
283     else:
284     return 0.
285    
286 gross 3478 def getBHP(self):
287     """
288     return bottom-hole pressure
289    
290     :note: needs to be overwritten
291     """
292     raise NotImplementedError
293    
294     def getPhase(self):
295     """
296     returns the pahse the well works on
297     """
298 gross 3484 return self.__phase
299    
300     def isOpen(self, t):
301     """
302     returns True is the well is open at time t
303     """
304     out = False
305     t0=min(t, self.__schedule[0])
306     i=0
307     while t < t0:
308     if out:
309     out=False
310     else:
311     out=True
312     i+=1
313     if i < len(self.__schedule[i]):
314     t0=self.__schedule[i]
315     else:
316     t0=t
317     return out
318 gross 3458
319 gross 3484 class VerticalPeacemanWell(Well):
320 gross 3478 """
321     defines a well using Peaceman's formula
322     """
323 gross 3484 def __init__(self,name, schedule = [0.], BHP=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
324     perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
325     """
326     set up well
327    
328     :param BHP: ottom-hole pressure
329     :param Q: flow rate
330     :param r: well radius
331     :param X: location
332     :param D: dimension of well block
333     :param perm: permeabilities
334     :param s: skin factor
335     """
336     Well.__init__(self, Q=Q, schedule=schedule, phase=phase,)
337     r_el=0.28 * ( (perm[1]/perm[0])**0.5 * D[0]**2 + (perm[0]/perm[1])**0.5 * D[1]**2 ) \
338     / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
339     self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
340     self.__r = r
341 gross 3458
342 gross 3484 self.__BHP = BHP
343     self.__D=D
344     self.r_el=r_el
345     self.X0=X0
346    
347     def getGeometry(self):
348     """
349     returns a function representing the geometry of the well.
350     typically a Gaussian profile around the location of the well.
351    
352     :note: needs to be overwritten
353     """
354     raise NotImplementedError
355    
356     def getProductivityIndex(self,t):
357     """
358     returns the productivity index of the well.
359     typically a Gaussian profile around the location of the well.
360     """
361     return self.__PI
362    
363     def getBHP(self):
364     """
365     return bottom-hole pressure
366     """
367     return self.__BHP
368    
369 gross 3478 class DualPorosity(object):
370     """
371     generic dual porosity model
372     """
373     def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
374     S_fg=None, S_mg=None,
375     perm_f_0=None, perm_f_1=None, perm_f_2=None,
376     perm_m_0=None, perm_m_1=None, perm_m_2=None,
377     k_w =None, k_g=None, mu_w =None, mu_g =None,
378     rho_w =None, rho_g=None,
379     wells=[]):
380     """
381     set up
382    
383     :param domain: domain
384     :type domain: `esys.escript.Domain`
385     :param phi_f: porosity of the fractured rock as function of the gas pressure
386     :type phi_f: `MaterialPropertyWithDifferential`
387     :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
388     :type phi_m: `MaterialPropertyWithDifferential` or None
389     :param phi: total porosity if None phi=1.
390     :type phi: scalar or None
391 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
392 gross 3478 :type S_fg: `MaterialPropertyWithDifferential`
393     :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
394     :type S_mg: `MaterialPropertyWithDifferential`
395     :param perm_f_0: permeability in the x_0 direction in the fractured media
396     :type perm_f_0: scalar
397     :param perm_f_1: permeability in the x_1 direction in the fractured media
398     :type perm_f_1: scalar
399     :param perm_f_2: permeability in the x_2 direction in the fractured media
400     :type perm_f_2: scalar
401     :param perm_m_0: permeability in the x_0 direction in the coal matrix
402     :type perm_m_0: scalar
403     :param perm_m_1: permeability in the x_1 direction in the coal matrix
404     :type perm_m_1: scalar
405     :param perm_m_2: permeability in the x_2 direction in the coal matrix
406     :type perm_m_2: scalar
407     :param k_w: relative permeability of water as function of water saturation
408     :type k_w: `MaterialProperty`
409     :param k_g: relative permeability of gas as function of gas saturation
410     :type k_g: `MaterialProperty`
411     :param mu_w: viscosity of water as function of water pressure
412     :type mu_w: `MaterialProperty`
413     :param mu_g: viscosity of gas as function of gas pressure
414     :type mu_g: `MaterialProperty`
415     :param rho_w: density of water as function of water pressure
416     :type rho_w: `MaterialPropertyWithDifferential`
417     :param rho_g: density of gas as function of gas pressure
418     :type rho_g: `MaterialPropertyWithDifferential`
419     :param wells : list of wells
420     :type wells: list of `Well`
421     """
422 gross 3458
423 gross 3478 self.domain = domain
424     self.phi_f = phi_f
425     self.phi_m = phi_m
426     self.phi = phi
427     self.S_fg = S_fg
428     self.S_mg = S_mg
429     self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
430     self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
431     self.k_w = k_w
432     self.k_g = k_g
433     self.mu_w = mu_w
434     self.mu_g = mu_g
435     self.rho_w = rho_w
436     self.rho_g = rho_g
437     self.wells=wells
438    
439 gross 3458
440    
441    
442    
443 gross 3478 class PorosityOneHalfModel(DualPorosity):
444     """
445     Model for gas and water in double prosity model tracking water and gas
446     pressure in the fractured rock and gas concentration in the coal matrix.
447     This corresponds to the coal bed methan model in the ECLIPSE code.
448     """
449 gross 3458
450 gross 3484 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
451 gross 3478 perm_f_0=None, perm_f_1=None, perm_f_2=None,
452     k_w =None, k_g=None, mu_w =None, mu_g =None,
453     rho_w =None, rho_g=None,
454     wells=[]):
455     """
456     set up
457    
458     :param domain: domain
459     :type domain: `esys.escript.Domain`
460     :param phi_f: porosity of the fractured rock as function of the gas pressure
461     :type phi_f: `MaterialPropertyWithDifferential`
462     :param phi: total porosity if None phi=1.
463     :type phi: scalar or None
464     :param L_g: gas adsorption as function of gas pressure
465     :type L_g: `MaterialPropertyWithDifferential`
466 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
467 gross 3478 :type S_fg: `MaterialPropertyWithDifferential`
468     :param perm_f_0: permeability in the x_0 direction in the fractured media
469     :type perm_f_0: scalar
470     :param perm_f_1: permeability in the x_1 direction in the fractured media
471     :type perm_f_1: scalar
472     :param perm_f_2: permeability in the x_2 direction in the fractured media
473     :type perm_f_2: scalar
474     :param k_w: relative permeability of water as function of water saturation
475     :type k_w: `MaterialProperty`
476     :param k_g: relative permeability of gas as function of gas saturation
477     :type k_g: `MaterialProperty`
478     :param mu_w: viscosity of water as function of water pressure
479     :type mu_w: `MaterialProperty`
480     :param mu_g: viscosity of gas as function of gas pressure
481     :type mu_g: `MaterialProperty`
482     :param rho_w: density of water as function of water pressure
483     :type rho_w: `MaterialPropertyWithDifferential`
484     :param rho_g: density of gas as function of gas pressure
485     :type rho_g: `MaterialPropertyWithDifferential`
486     :param wells : list of wells
487     :type wells: list of `Well`
488     """
489    
490     DualPorosity.__init__(self, domain,
491     phi_f=phi_f, phi_m=None, phi=phi,
492 gross 3484 S_fg=None, S_mg=None,
493 gross 3478 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
494     perm_m_0=None , perm_m_1=None, perm_m_2=None,
495     k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
496     rho_w =rho_w, rho_g=rho_g,
497     wells=wells)
498     self.L_g=L_g
499     self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
500    
501 gross 3458
502 gross 3478 def getPDEOptions(self):
503     """
504     returns the `SolverOptions` of the underlying PDE
505     """
506     return self.__pde.getSolverOptions()
507    
508 gross 3484 def setInitialState(p=1.*U.atm, S_fg=0, C_mg=None):
509     """
510     sets the initial state
511    
512     :param p: pressure
513     :param S_fg: gas saturation in fractured rock
514     :param C_mg: gas concentration in coal matrix. if not given it is calculated
515     using the gas adsorption curve.
516     """
517    
518     self.u=self.__pde.getNewCoefficient("u")
519     u[0]=p
520     u[1]=S_fg
521     if C_mg == None:
522     u[2]= self.L_g(p)
523     else:
524     u[2]=C_mg
525    
526 gross 3478 def solvePDE(self):
527    
528 gross 3484 p_f=self.u[0]
529     S_fg=self.u[1]
530 gross 3478 C_mg=self.u[3]
531    
532    
533 gross 3484 p_f_old=self.u_old[0]
534     S_fg_old=self.u_old[1]
535 gross 3478 C_mg_old=self.u_old[3]
536 gross 3458
537    
538 gross 3484 phi_f =self.phi_f.getValue(S_fg)
539     dphi_fdp=self.phi_f.getValueDifferential(S_fg)
540 gross 3478
541     S_fw= 1-S_fg
542    
543 gross 3484 rho_fw = self.rho_w.getValue(p_f)
544     drho_fwdp = self.rho_w.getValueDifferential( p_f)
545     rho_fg = self.rho_g.getValue(p_f)
546     drho_fgdp = self.rho_g.getValueDifferential(p_f)
547 gross 3478
548 gross 3484 L_g = self.getValue(p_f)
549     dL_gdp = self.rho_w.getValueDifferential(p_f)
550 gross 3478
551     A_fw = rho_fw * k_fw/mu_fw
552     A_fg = rho_fg * k_fg/mu_fg
553    
554     m = whereNegative(L_g - C_mg)
555     B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )
556    
557    
558 gross 3484 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
559     E_fps= - phi_f * rho_fw
560     E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f * drho_fgdp )
561     E_fss= phi_f * rho_fg
562 gross 3458
563    
564 gross 3478
565     Q_w=0.
566     Q_g=0.
567     R_w=0.
568     R_g=0.
569    
570     for I in self.wells:
571     chi_I= I.getGeometry()
572     Pi_I = I.getProductivityIndex()
573     P_I = I.getBHP()
574     Q_I = I.getFlowRate()
575    
576     if self.getPhase == Well.WATER:
577     Q_w = Q_w + Q_I * chi_I
578     Q_g = Q_g + Pi_I * P_I * chi_I
579     R_g = R_g + Pi_I * chi_I
580     else:
581     Q_g = Q_g + Q_I * chi_I
582     Q_w = Q_w + Pi_I * P_I * chi_I
583     R_w = R_w + Pi_I * chi_I
584    
585 gross 3484 F_fp_Y = A_fw * Q_w
586     F_fs_Y = A_fg * Q_g
587     D_fss = - A_fg * R_g
588     D_fpp = - A_fw * R_w
589 gross 3478
590     if self.domain.getDim() == 3:
591 gross 3484 F_fp_X = A_fw * self.g * rho_fw * self.perm_f[2] * kronecker(self.domain)[2]
592     F_fs_X = A_fg * self.g * rho_fg * self.perm_f[2] * kronecker(self.domain)[2]
593 gross 3478 else:
594 gross 3484 F_fp_X = 0 * kronecker(self.domain)[1]
595     F_fs_X = 0 * kronecker(self.domain)[1]
596 gross 3478
597 gross 3484 F_mg = B * ( L_g - dL_gdp * S_fg )
598 gross 3458
599 gross 3478
600     D=self.__pde.getNewCoefficient("D")
601     A=self.__pde.getNewCoefficient("A")
602     Y=self.__pde.getNewCoefficient("Y")
603     X=self.__pde.getNewCoefficient("X")
604    
605     dtXI = dt*self.XI
606     dtcXI = dt*(1-self.XI)
607    
608 gross 3484 D[0,0]=E_fpp + dtXI * D_fpp
609     D[0,1]=E_fps
610     D[1,0]=E_fsp
611     D[1,1]=E_fss + dtXI * D_fss
612 gross 3478 D[1,2]=rho_fg
613 gross 3484 D[2,1]= dtXI * B * dL_gdp
614 gross 3478 D[2,2]= 1 + dtXI * B
615    
616 gross 3484 H_fw = dtcXI * A_fw * grad( p_f_old)
617     H_fg = dtcXI * A_fg * grad(S_fg_old)
618 gross 3478
619     for i in range(self.domain.getDim()):
620     A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
621     A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
622    
623 gross 3484 X[0,i]= dt * F_fp_X - H_fw * self.perm_f[i]
624     X[1,i]= dt * F_fs_X - H_fg * self.perm_f[i]
625 gross 3478
626 gross 3484 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old + dt * F_fp_Y - dtcXI * D_fpp * p_f_old
627     Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + rho_fg * C_mg_old + dt * F_fs_Y - dtcXI * D_fss * S_fg_old
628     Y[2] = C_mg_old + dt * F_mg_Y - dtcXI * ( dL_gdp * S_fg_old + C_mg_old)
629 gross 3478
630     self.__pde.setValue(A=A, D=D, X=X, Y=Y)
631    
632     u2 = self.__pde.getSolution()
633     return u2

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26