/[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 3485 - (show annotations)
Thu Mar 24 22:44:40 2011 UTC (8 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 19462 byte(s)
some compile fixes
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 self.u=self.__pde.getNewCoefficient("u")
518 u[0]=p
519 u[1]=S_fg
520 if C_mg == None:
521 u[2]= self.L_g(p)
522 else:
523 u[2]=C_mg
524
525 def solvePDE(self):
526
527 p_f=self.u[0]
528 S_fg=self.u[1]
529 C_mg=self.u[3]
530
531
532 p_f_old=self.u_old[0]
533 S_fg_old=self.u_old[1]
534 C_mg_old=self.u_old[3]
535
536
537 phi_f =self.phi_f.getValue(S_fg)
538 dphi_fdp=self.phi_f.getValueDifferential(S_fg)
539
540 S_fw= 1-S_fg
541
542 rho_fw = self.rho_w.getValue(p_f)
543 drho_fwdp = self.rho_w.getValueDifferential( p_f)
544 rho_fg = self.rho_g.getValue(p_f)
545 drho_fgdp = self.rho_g.getValueDifferential(p_f)
546
547 L_g = self.getValue(p_f)
548 dL_gdp = self.rho_w.getValueDifferential(p_f)
549
550 A_fw = rho_fw * k_fw/mu_fw
551 A_fg = rho_fg * k_fg/mu_fg
552
553 m = whereNegative(L_g - C_mg)
554 B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )
555
556
557 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
558 E_fps= - phi_f * rho_fw
559 E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f * drho_fgdp )
560 E_fss= phi_f * rho_fg
561
562
563
564 Q_w=0.
565 Q_g=0.
566 R_w=0.
567 R_g=0.
568
569 for I in self.wells:
570 chi_I= I.getGeometry()
571 Pi_I = I.getProductivityIndex()
572 P_I = I.getBHP()
573 Q_I = I.getFlowRate()
574
575 if self.getPhase == Well.WATER:
576 Q_w = Q_w + Q_I * chi_I
577 Q_g = Q_g + Pi_I * P_I * chi_I
578 R_g = R_g + Pi_I * chi_I
579 else:
580 Q_g = Q_g + Q_I * chi_I
581 Q_w = Q_w + Pi_I * P_I * chi_I
582 R_w = R_w + Pi_I * chi_I
583
584 F_fp_Y = A_fw * Q_w
585 F_fs_Y = A_fg * Q_g
586 D_fss = - A_fg * R_g
587 D_fpp = - A_fw * R_w
588
589 if self.domain.getDim() == 3:
590 F_fp_X = A_fw * self.g * rho_fw * self.perm_f[2] * kronecker(self.domain)[2]
591 F_fs_X = A_fg * self.g * rho_fg * self.perm_f[2] * kronecker(self.domain)[2]
592 else:
593 F_fp_X = 0 * kronecker(self.domain)[1]
594 F_fs_X = 0 * kronecker(self.domain)[1]
595
596 F_mg = B * ( L_g - dL_gdp * S_fg )
597
598
599 D=self.__pde.getNewCoefficient("D")
600 A=self.__pde.getNewCoefficient("A")
601 Y=self.__pde.getNewCoefficient("Y")
602 X=self.__pde.getNewCoefficient("X")
603
604 dtXI = dt*self.XI
605 dtcXI = dt*(1-self.XI)
606
607 D[0,0]=E_fpp + dtXI * D_fpp
608 D[0,1]=E_fps
609 D[1,0]=E_fsp
610 D[1,1]=E_fss + dtXI * D_fss
611 D[1,2]=rho_fg
612 D[2,1]= dtXI * B * dL_gdp
613 D[2,2]= 1 + dtXI * B
614
615 H_fw = dtcXI * A_fw * grad( p_f_old)
616 H_fg = dtcXI * A_fg * grad(S_fg_old)
617
618 for i in range(self.domain.getDim()):
619 A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
620 A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
621
622 X[0,i]= dt * F_fp_X - H_fw * self.perm_f[i]
623 X[1,i]= dt * F_fs_X - H_fg * self.perm_f[i]
624
625 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old + dt * F_fp_Y - dtcXI * D_fpp * p_f_old
626 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
627 Y[2] = C_mg_old + dt * F_mg_Y - dtcXI * ( dL_gdp * S_fg_old + C_mg_old)
628
629 self.__pde.setValue(A=A, D=D, X=X, Y=Y)
630
631 u2 = self.__pde.getSolution()
632 return u2

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26