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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3484 - (show annotations)
Thu Mar 24 05:42:07 2011 UTC (8 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 19464 byte(s)
some work on coal gas
1 #######################################################
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 from esys.escript.linearPDEs import LinearPDE
24 from esys.escript import unitsSI as U
25 from esys.escript import sqrt, log
26 import math
27
28 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
53 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
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
84 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 """
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 MaterialPropertyWithDifferential.__init__(self)
184 if len(x) != len(y):
185 raise ValueError,"length of interpolation nodes and value lists need to be identical."
186 if len(x) < 1 :
187 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
198 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
221 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 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 def __init__(self, Q=0., schedule=[0.], phase="Water", *args, **kwargs):
252 """
253 set up well
254 """
255 self.__schedule=schedule
256 self.__phase=phase
257 self.__Q = Q
258
259 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
268 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 def getFlowRate(self, t):
278 """
279 returns the flow rate
280 """
281 if self.isOpen(t):
282 return self.__Q
283 else:
284 return 0.
285
286 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 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
319 class VerticalPeacemanWell(Well):
320 """
321 defines a well using Peaceman's formula
322 """
323 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
342 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 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 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
392 :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
423 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
440
441
442
443 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
450 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
451 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 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
467 :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 S_fg=None, S_mg=None,
493 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
502 def getPDEOptions(self):
503 """
504 returns the `SolverOptions` of the underlying PDE
505 """
506 return self.__pde.getSolverOptions()
507
508 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 def solvePDE(self):
527
528 p_f=self.u[0]
529 S_fg=self.u[1]
530 C_mg=self.u[3]
531
532
533 p_f_old=self.u_old[0]
534 S_fg_old=self.u_old[1]
535 C_mg_old=self.u_old[3]
536
537
538 phi_f =self.phi_f.getValue(S_fg)
539 dphi_fdp=self.phi_f.getValueDifferential(S_fg)
540
541 S_fw= 1-S_fg
542
543 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
548 L_g = self.getValue(p_f)
549 dL_gdp = self.rho_w.getValueDifferential(p_f)
550
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 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
563
564
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 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
590 if self.domain.getDim() == 3:
591 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 else:
594 F_fp_X = 0 * kronecker(self.domain)[1]
595 F_fs_X = 0 * kronecker(self.domain)[1]
596
597 F_mg = B * ( L_g - dL_gdp * S_fg )
598
599
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 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 D[1,2]=rho_fg
613 D[2,1]= dtXI * B * dL_gdp
614 D[2,2]= 1 + dtXI * B
615
616 H_fw = dtcXI * A_fw * grad( p_f_old)
617 H_fg = dtcXI * A_fg * grad(S_fg_old)
618
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 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
626 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
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