/[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 3497 - (show annotations)
Wed Apr 6 06:11:56 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 26278 byte(s)

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.pdetools import Locator
26 from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip
27 from esys.weipa import saveVTK
28 import math
29
30 class MaterialProperty(object):
31 """
32 generic class for material properties depending on one (or several) status
33 variable(s)
34 """
35 def __init__(self, *args, **kwargs):
36 """
37 set up
38 """
39 pass
40 def __call__(self, *args, **kwargs):
41 """
42 return value of material property for given status
43 """
44 return self.getValue( *args, **kwargs)
45
46
47 def getValue(self, *args, **kwargs):
48 """
49 return value of material property for given status
50
51 :remark: needs to be overwritten
52 """
53 raise NotImplementedError
54
55 class MaterialPropertyWithDifferential(MaterialProperty):
56 """
57 generic class for material properties depending on one (or several) status
58 variable(s) where the derivative of the property with respect to the
59 status variables is available
60 """
61 def getValueDifferential(self, *args, **kwargs):
62 """
63 returns the value of the derivative of material property for given status variable
64
65 :remark: needs to be overwritten
66 """
67 raise NotImplementedError
68
69 class Porosity(MaterialPropertyWithDifferential):
70 """
71 defines porosity phi as function of pressure
72
73 phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
74
75 phi_p_ref is claculted from the initial porosity phi_0 at pressure p_0
76 """
77 def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
78 """
79 """
80 self.phi_p_ref=1 # will be overwritten later
81 self.p_ref=p_ref
82 self.C=C
83 # reset phi_p_ref to get phi(p_0)=phi_0
84 self.phi_p_ref=phi_0/self.getValue(p_0)
85
86 def getValue(self, p):
87 """
88 returns the porosity for given pressure p
89 """
90 X= self.C * ( p - self.p_ref )
91 return self.phi_p_ref * (1. + X * (1. + X/2 ) )
92
93 def getValueDifferential(self, p):
94 """
95 returns the porosity for given pressure p
96 """
97 X= self.C * ( p - self.p_ref )
98 return self.phi_p_ref * self.C * (1. + X )
99
100
101 class WaterDensity(MaterialPropertyWithDifferential):
102 """
103 set water density as
104
105 rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
106
107 with rho_ref =rho_s/B_ref * gravity
108 """
109 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):
110 self.rho_0 = rho_s * gravity
111 self.rho_ref = self.rho_0/B_ref
112 self.C=C
113 self.p_ref=p_ref
114
115 def getValue(self, p):
116 """
117 returns the density for given pressure p
118 """
119 X= self.C * ( p - self.p_ref )
120 return self.rho_ref * (1+ X * (1+ X/2) )
121
122 def getValueDifferential(self, p):
123 """
124 """
125 X= self.C * ( p - self.p_ref )
126 return self.rho_ref * self.C * (1+ X)
127
128 def getFormationFactor(self, p):
129 return self.rho_0/self.getValue(p)
130
131
132
133 class WaterViscosity(MaterialProperty):
134 """
135 defines viscosity of water
136
137 mu=mu_ref /(1 + X + X*X/2)
138
139 with X=-C*(p-p_ref)
140 """
141
142 def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
143 """
144 set up
145 """
146 self.mu_ref = mu_ref
147 self.p_ref = p_ref
148 self.C=C
149 def getValue(self, p):
150 """
151 returns the viscosity for given pressure p
152 """
153 X= -self.C * ( p - self.p_ref )
154 return self.mu_ref/(1+ X * (1+ X/2) )
155
156 class GasDensity(MaterialPropertyWithDifferential):
157 """
158 set water density as
159
160 rho = gravity * rho_air_s /B(p)
161
162 where B is given by an interpolation table
163 """
164 def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
165 self.rho_ref =rho_air_s * gravity
166 self.rho_0 =rho_air_s * gravity
167 self.tab=InterpolationTable(x=p, y=B)
168
169 def getValue(self, p):
170 """
171 returns the density for given pressure p
172 """
173 return self.rho_ref/self.getFormationFactor(p)
174
175 def getValueDifferential(self, p):
176 """
177 """
178 B = self.getFormationFactor(p)
179 dBdp = self.getFormationFactorDifferential(p)
180 return -self.rho_ref * dBdp/(B * B)
181
182 def getFormationFactor(self, p):
183 return self.tab.getValue(p)
184
185 def getFormationFactorDifferential(self, p):
186 return self.tab.getValueDifferential(p)
187
188
189 class InterpolationTable(MaterialPropertyWithDifferential):
190 """
191 a simple 1D interpolation table for escript Data object with non-equidistant nodes
192 """
193 def __init__(self,x=[], y=[], obeyBounds=True):
194 """
195 set up interpolation table. obeyBounds is set an exception is thrown if
196 the interpolation argument is below min(x) or larger than max(x). Otherwise
197 the value for x is set to y[0] for
198 """
199 MaterialPropertyWithDifferential.__init__(self)
200 if len(x) != len(y):
201 raise ValueError,"length of interpolation nodes and value lists need to be identical."
202 if len(x) < 1 :
203 raise ValueError,"length of interpolation nodes a list needs to at least one."
204
205 x_ref=x[0]
206 for i in range(1,len(x)):
207 if x_ref >= x[i]:
208 raise ValueError,"interpolation nodes need to be given in increasing order."
209 x_ref=x[i]
210 self.__x = x
211 self.__y = y
212 self.__obeyBounds = obeyBounds
213
214 def getValue(self, x):
215 """
216 returns the interpolated values of x
217 """
218 X=self.__x
219 Y=self.__y
220
221 x0=X[0]
222 m0=whereNegative(x-x0)
223 if self.__obeyBounds:
224 if sup(m0) > 0:
225 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
226 out=self.__x[0]
227 for i in range(1,len(X)):
228 z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
229 out = out * m0 + z * (1-m0)
230 m0=whereNonPositive(x-X[i])
231
232 if self.__obeyBounds:
233 if inf(m0) < 1 :
234 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
235 else:
236 out = out * m0 + V[-1] * (1-m0)
237 return out
238
239 def getValueDifferential(self, x):
240 X=self.__x
241 Y=self.__y
242
243 x0=X[0]
244 m0=whereNegative(x-x0)
245 if self.__obeyBounds:
246 if sup(m0) > 0:
247 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
248 out=0.
249 for i in range(1,len(X)):
250 z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
251 out = out * m0 + z * (1-m0)
252 m0=whereNonPositive(x-X[i])
253
254 if self.__obeyBounds:
255 if inf(m0) < 1:
256 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
257 else:
258 out = out * m0
259 return out
260
261
262 class Well(object):
263 """
264 generic well
265
266 :var WATER: phase identifier for water well
267 :var GAS: phase identifier for gas well
268 """
269 WATER="Water"
270 GAS="Gas"
271 def __init__(self, name, domain, Q=0., schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
272 """
273 set up well
274 """
275 self.__schedule=schedule
276 self.__phase=phase
277 self.__Q = Q
278 self.__BHP_limit=BHP_limit
279 self.name=name
280 self.domain=domain
281
282 def getGeometry(self):
283 """
284 returns a function representing the geometry of the well.
285 typically a Gaussian profile around the location of the well.
286
287 :note: needs to be overwritten
288 """
289 raise NotImplementedError
290
291 def getProductivityIndex(self):
292 """
293 returns the productivity index of the well.
294 typically a Gaussian profile around the location of the well.
295
296 :note: needs to be overwritten
297 """
298 raise NotImplementedError
299
300 def getFlowRate(self):
301 """
302 returns the flow rate
303 """
304 return self.__Q
305
306 def getBHPLimit(self):
307 """
308 return bottom-hole pressure limit
309
310 :note: needs to be overwritten
311 """
312 return self.__BHP_limit
313
314 def getBHP(self):
315 """
316 return bottom-hole pressure
317
318 :note: needs to be overwritten
319 """
320 return self.__BHP
321
322 def setBHP(self, BHP):
323 """
324 sets bottom-hole pressure
325 """
326 self.__BHP= BHP
327
328 def getPhase(self):
329 """
330 returns the pahse the well works on
331 """
332 return self.__phase
333
334 def isOpen(self, t):
335 """
336 returns True is the well is open at time t
337 """
338 out = False
339 t0=min(t, self.__schedule[0])
340 i=0
341 while t > t0:
342 if out:
343 out=False
344 else:
345 out=True
346 i+=1
347 if i < len(self.__schedule):
348 t0=self.__schedule[i]
349 else:
350 t0=t
351 return out
352
353 def setWaterProduction(self, v):
354 self.water_production=v
355 def getWaterProduction(self):
356 return self.water_production
357 def setGasProduction(self, v):
358 self.gas_production=v
359 def getGasProduction(self):
360 return self.gas_production
361
362 class VerticalPeacemanWell(Well):
363 """
364 defines a well using Peaceman's formula
365 """
366 def __init__(self,name, domain, schedule = [0.], BHP_limit=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
367 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
368 """
369 set up well
370
371 :param BHP: ottom-hole pressure
372 :param Q: flow rate
373 :param r: well radius
374 :param X: location
375 :param D: dimension of well block
376 :param perm: permeabilities
377 :param s: skin factor
378 """
379 Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
380 r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 + sqrt(perm[0]/perm[1]) * D[1]**2 )\
381 / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
382 self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
383 self.__r = r
384
385 self.__D=D
386 self.r_el=r_el
387 self.X0=X0[:self.domain.getDim()]
388
389 x=Function(domain).getX()
390 self.chi = 1.
391 for i in range(domain.getDim()):
392 self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)
393
394 self.chi*=1./(D[0]*D[1]*D[2])
395
396
397 #self.chi=0.00000001
398 def getLocation(self):
399 return self.X0
400
401 def getGeometry(self):
402 """
403 returns a function representing the geometry of the well.
404 typically a Gaussian profile around the location of the well.
405
406 :note: needs to be overwritten
407 """
408 return self.chi
409
410 def getProductivityIndex(self):
411 """
412 returns the productivity index of the well.
413 """
414 return self.__PI
415
416
417 class DualPorosity(object):
418 """
419 generic dual porosity model
420 """
421 def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
422 S_fg=None, S_mg=None,
423 perm_f_0=None, perm_f_1=None, perm_f_2=None,
424 perm_m_0=None, perm_m_1=None, perm_m_2=None,
425 k_w =None, k_g=None, mu_w =None, mu_g =None,
426 rho_w =None, rho_g=None,
427 wells=[]):
428 """
429 set up
430
431 :param domain: domain
432 :type domain: `esys.escript.Domain`
433 :param phi_f: porosity of the fractured rock as function of the gas pressure
434 :type phi_f: `MaterialPropertyWithDifferential`
435 :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
436 :type phi_m: `MaterialPropertyWithDifferential` or None
437 :param phi: total porosity if None phi=1.
438 :type phi: scalar or None
439 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
440 :type S_fg: `MaterialPropertyWithDifferential`
441 :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
442 :type S_mg: `MaterialPropertyWithDifferential`
443 :param perm_f_0: permeability in the x_0 direction in the fractured media
444 :type perm_f_0: scalar
445 :param perm_f_1: permeability in the x_1 direction in the fractured media
446 :type perm_f_1: scalar
447 :param perm_f_2: permeability in the x_2 direction in the fractured media
448 :type perm_f_2: scalar
449 :param perm_m_0: permeability in the x_0 direction in the coal matrix
450 :type perm_m_0: scalar
451 :param perm_m_1: permeability in the x_1 direction in the coal matrix
452 :type perm_m_1: scalar
453 :param perm_m_2: permeability in the x_2 direction in the coal matrix
454 :type perm_m_2: scalar
455 :param k_w: relative permeability of water as function of water saturation
456 :type k_w: `MaterialProperty`
457 :param k_g: relative permeability of gas as function of gas saturation
458 :type k_g: `MaterialProperty`
459 :param mu_w: viscosity of water as function of water pressure
460 :type mu_w: `MaterialProperty`
461 :param mu_g: viscosity of gas as function of gas pressure
462 :type mu_g: `MaterialProperty`
463 :param rho_w: density of water as function of water pressure
464 :type rho_w: `MaterialPropertyWithDifferential`
465 :param rho_g: density of gas as function of gas pressure
466 :type rho_g: `MaterialPropertyWithDifferential`
467 :param wells : list of wells
468 :type wells: list of `Well`
469 """
470
471 self.domain = domain
472 self.phi_f = phi_f
473 self.phi_m = phi_m
474 self.phi = phi
475 self.S_fg = S_fg
476 self.S_mg = S_mg
477 self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
478 self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
479 self.k_w = k_w
480 self.k_g = k_g
481 self.mu_w = mu_w
482 self.mu_g = mu_g
483 self.rho_w = rho_w
484 self.rho_g = rho_g
485 self.wells=wells
486 self.t =0
487
488 self.__iter_max=1
489 self.__rtol=1.e-4
490 self.verbose=False
491 self.XI=0.5
492 def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
493 """
494 sets parameters to control iteration process
495 """
496 if iter_max !=None: self.__iter_max=iter_max
497 if rtol !=None: self.__rtol = rtol
498 if verbose !=None: self.verbose=verbose
499 if xi !=None: self.XI=xi
500
501 def update(self, dt):
502 self.u, self.u_old = self.u.copy(), self.u
503 n=0
504 rerr=1.
505 while n < self.__iter_max and rerr > self.__rtol:
506 u=self.solvePDE(dt)
507 norm_u=Lsup(u)
508 norm_e=Lsup(u-self.u)
509
510 if norm_u > 0:
511 rerr=norm_e/norm_u
512 else:
513 rerr=norm_e
514 if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)
515 n+=1
516 self.u=u
517 print "iteration completed."
518 self.t+=dt
519
520
521 class PorosityOneHalfModel(DualPorosity):
522 """
523 Model for gas and water in double prosity model tracking water and gas
524 pressure in the fractured rock and gas concentration in the coal matrix.
525 This corresponds to the coal bed methan model in the ECLIPSE code.
526 """
527
528 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
529 perm_f_0=None, perm_f_1=None, perm_f_2=None,
530 k_w =None, k_g=None, mu_w =None, mu_g =None,
531 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
532 wells=[]):
533 """
534 set up
535
536 :param domain: domain
537 :type domain: `esys.escript.Domain`
538 :param phi_f: porosity of the fractured rock as function of the gas pressure
539 :type phi_f: `MaterialPropertyWithDifferential`
540 :param phi: total porosity if None phi=1.
541 :type phi: scalar or None
542 :param L_g: gas adsorption as function of gas pressure
543 :type L_g: `MaterialPropertyWithDifferential`
544 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
545 :type S_fg: `MaterialPropertyWithDifferential`
546 :param perm_f_0: permeability in the x_0 direction in the fractured media
547 :type perm_f_0: scalar
548 :param perm_f_1: permeability in the x_1 direction in the fractured media
549 :type perm_f_1: scalar
550 :param perm_f_2: permeability in the x_2 direction in the fractured media
551 :type perm_f_2: scalar
552 :param k_w: relative permeability of water as function of water saturation
553 :type k_w: `MaterialProperty`
554 :param k_g: relative permeability of gas as function of gas saturation
555 :type k_g: `MaterialProperty`
556 :param mu_w: viscosity of water as function of water pressure
557 :type mu_w: `MaterialProperty`
558 :param mu_g: viscosity of gas as function of gas pressure
559 :type mu_g: `MaterialProperty`
560 :param rho_w: density of water as function of water pressure
561 :type rho_w: `MaterialPropertyWithDifferential`
562 :param rho_g: density of gas as function of gas pressure
563 :type rho_g: `MaterialPropertyWithDifferential`
564 :param wells : list of wells
565 :type wells: list of `Well`
566 :param sigma: shape factor for gas matrix diffusion
567 :param A_mg: diffusion constant for gas adsorption
568 :param f_rg: gas re-adsorption factor
569 """
570
571 DualPorosity.__init__(self, domain,
572 phi_f=phi_f, phi_m=None, phi=phi,
573 S_fg=None, S_mg=None,
574 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
575 perm_m_0=None , perm_m_1=None, perm_m_2=None,
576 k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
577 rho_w =rho_w, rho_g=rho_g,
578 wells=wells)
579 self.L_g=L_g
580 self.sigma = sigma
581 self.A_mg = A_mg
582 self.f_rg = f_rg
583 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
584
585
586 def getPDEOptions(self):
587 """
588 returns the `SolverOptions` of the underlying PDE
589 """
590 return self.__pde.getSolverOptions()
591
592 def getState(self):
593 return self.u[0], self.u[1], self.u[2]
594
595 def getOldState(self):
596 return self.u_old[0], self.u_old[1], self.u_old[2]
597
598 def setInitialState(self, p=1.*U.atm, S_fg=0, C_mg=None):
599 """
600 sets the initial state
601
602 :param p: pressure
603 :param S_fg: gas saturation in fractured rock
604 :param C_mg: gas concentration in coal matrix. if not given it is calculated
605 using the gas adsorption curve.
606 """
607 self.u=self.__pde.createSolution()
608 self.u[0]=p
609 self.u[1]=S_fg
610 if C_mg == None:
611 self.u[2]= self.L_g(p)*self.rho_g.getFormationFactor(p)
612 else:
613 self.u[2]=C_mg
614
615 def solvePDE(self, dt):
616
617 C_couple=1
618
619 p_f, S_fg, C_mg=self.getState()
620 p_f_old, S_fg_old, C_mg_old=self.getOldState()
621
622 S_fw=1-S_fg
623
624 if self.verbose:
625 print "p_f range = ",inf(p_f),sup(p_f)
626 print "S_fg range = ",inf(S_fg),sup(S_fg)
627 print "S_fw range = ",inf(S_fw),sup(S_fw)
628 print "C_mg range = ",inf(C_mg),sup(C_mg)
629
630 k_fw=self.k_w(S_fw)
631 if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
632
633 k_fg=self.k_g(S_fg)
634 if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
635
636 mu_fw=self.mu_w(p_f)
637 if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
638
639 mu_fg=self.mu_g(p_f)
640 if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
641
642
643 phi_f =self.phi_f.getValue(p_f)
644 dphi_fdp=self.phi_f.getValueDifferential(p_f)
645 if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
646
647 rho_fw = self.rho_w.getValue(p_f)
648 drho_fwdp = self.rho_w.getValueDifferential(p_f)
649 if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
650
651 rho_fg = self.rho_g.getValue(p_f)
652 drho_fgdp = self.rho_g.getValueDifferential(p_f)
653 if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
654
655 L_g_0 = self.L_g(p_f)
656 dL_g_0dp = self.L_g.getValueDifferential(p_f)
657 FF_g=self.rho_g.getFormationFactor(p_f)
658 dFF_gdp=self.rho_g.getFormationFactorDifferential(p_f)
659
660 L_g = L_g_0 * FF_g
661 dL_gdp = 0 * (dL_g_0dp * FF_g + L_g_0 * dFF_gdp)
662
663 if self.verbose:
664 print "L_g_0 range = ",inf(L_g_0),sup(L_g_0)
665 print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
666
667
668 A_fw = rho_fw * k_fw/mu_fw
669 A_fg = rho_fg * k_fg/mu_fg
670
671 m = whereNegative(L_g - C_mg)
672 B = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
673 print "XXXX B =",B
674 print "XXXX self.sigma =",self.sigma
675 print "XXXX self.A_mg =",self.A_mg
676
677 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
678 E_fps= - phi_f * rho_fw
679 E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f * drho_fgdp )
680 E_fss= phi_f * rho_fg
681
682 F_fw=0.
683 F_fg=0.
684 D_fw=0.
685 D_fg=0.
686
687 for I in self.wells:
688 chi_I= I.getGeometry()
689 loc=Locator(Function(self.domain),I.getLocation())
690 Pi_I = I.getProductivityIndex()
691 A_fw_I= loc(A_fw)
692 A_fg_I= loc(A_fg)
693 BHP_limit_I=I.getBHPLimit()
694
695 if I.isOpen(self.t+dt):
696 if self.verbose: print "well %s is open."%I.name
697 if I.getPhase() == Well.WATER:
698 q_I = self.rho_w.rho_0 * I.getFlowRate()
699 p_I_ref=loc(p_f)-q_I/(A_fw_I * Pi_I)
700 else:
701 q_I = self.rho_g.rho_0 * I.getFlowRate()
702 p_I_ref=loc(p_f)-q_I/(A_fg_I * Pi_I)
703
704 print "ZZZ =",loc(p_f), q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)
705
706 if BHP_limit_I > p_I_ref:
707 BHP_I=BHP_limit_I
708 D_fw = D_fw + Pi_I * A_fw_I * chi_I
709 F_fw = F_fw - Pi_I * A_fw_I * BHP_limit_I *chi_I
710 D_fg = D_fg + Pi_I * A_fg_I * chi_I
711 F_fg = F_fg - Pi_I * A_fg_I * BHP_limit_I *chi_I
712 else:
713 if I.getPhase() == Well.WATER:
714 F_fw = F_fw + q_I * chi_I
715 F_fg = F_fg + A_fg_I/A_fw_I * q_I *chi_I
716 else:
717 F_fg = F_fg + q_I * chi_I
718 F_fw = F_fw + A_fw_I/A_fg_I * q_I *chi_I
719 BHP_I=p_I_ref
720 else:
721 if self.verbose: print "well %s is closed."%I.name
722 BHP_I=loc(p_f)
723
724 if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
725 I.setBHP(BHP_I)
726
727 F_fp_Y = - F_fw
728 F_fs_Y = - F_fg
729 D_fpp = D_fw
730 D_fsp = D_fg
731
732
733 if self.domain.getDim() == 3:
734 F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
735 F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
736 else:
737 F_fp_X = 0 * kronecker(self.domain)[1]
738 F_fs_X = 0 * kronecker(self.domain)[1]
739
740 #F_mg_Y = B * ( L_g - dL_gdp * p_f )
741 F_mg_Y = B * L_g
742
743
744 D=self.__pde.createCoefficient("D")
745 A=self.__pde.createCoefficient("A")
746 Y=self.__pde.createCoefficient("Y")
747 X=self.__pde.createCoefficient("X")
748
749 dtXI = dt*self.XI
750 dtcXI = dt*(1-self.XI)
751
752 D[0,0]=E_fpp + dtXI * D_fpp
753 D[0,1]=E_fps
754 D[1,0]=E_fsp + dtXI * D_fsp
755 D[1,1]=E_fss
756 D[1,2]=rho_fg * C_couple
757 #D[2,0]= - dtXI * B * dL_gdp
758 D[2,2]= 1 + dtXI * B
759
760
761 for i in range(self.domain.getDim()):
762 A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
763 A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
764
765 g= grad(p_f_old) * dtcXI * self.perm_f[0:self.domain.getDim()]
766 X[0,:]= - A_fw * g + dt * F_fp_X
767 X[1,:]= - A_fg * g + dt * F_fs_X
768
769 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old + dt * F_fp_Y - dtcXI * D_fpp * p_f_old
770 Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + C_couple * rho_fg * C_mg_old + dt * F_fs_Y - dtcXI * D_fsp * p_f_old
771 Y[2] = C_mg_old + dt * F_mg_Y - dtcXI * B * C_mg_old # - dL_gdp * p_f_old)
772
773 print "HHH Y[0] =", Y[0]
774 print "HHH A_fw = ",A_fw
775 print "HHH A_fg= ",A_fg
776 print "HHH F_fg = ",F_fg
777 print "HHH F_fw = ",F_fw
778 print "HHH perm_f = ",self.perm_f
779 print "HHH k = ",(A_fw*self.perm_f[0])/D[0,0]
780 print "HHH k = ",(A_fw*self.perm_f[1])/D[0,0]
781 print "HHH X[0,:]= ",X[0,:]
782 print "HHH D[0,0] = ",D[0,0]
783 print "HHH D[0,1] = ",D[0,1]
784 print "HHH D[0,2] = ",D[0,2]
785
786 print "hhh B= ",B
787
788
789
790
791 #print "HHH Y[1] = ",Y[1]
792 #print "HHH X[1,:] = ",X[1,:]
793 #print "HHH D[1,0] = ",D[1,0]
794 #print "HHH D[1,1] = ",D[1,1]
795 #print "HHH D[1,2] = ",D[1,2]
796 #print "HHH A_fg =", A_fg
797
798
799
800
801 #1/0
802 self.__pde.setValue(A=A, D=D, X=X, Y=Y)
803
804 u2 = self.__pde.getSolution()
805 # this is deals with round-off errors to maintain physical meaningful values
806 # we should do this in a better way to detect values that are totally wrong
807 u2[1]=clip(u2[1],minval=0, maxval=1)
808 u2[2]=clip(u2[2],minval=0)
809
810 # update well production
811 p_f=u2[0]
812 for I in self.wells:
813 loc=Locator(Function(self.domain),I.getLocation())
814 Pi_I = I.getProductivityIndex()
815 q_I = I.getFlowRate()
816 A_fw_I= loc(A_fw)
817 A_fg_I= loc(A_fg)
818 p_f_I=loc(A_fg)
819 BHP_limit_I=I.getBHPLimit()
820 BHP=I.getBHP()
821
822
823 rho_fw_I = loc(rho_fw)
824 rho_fg_I = loc(rho_fg)
825 A_fg = rho_fg * k_fg/mu_fg
826
827 if I.isOpen(self.t+dt):
828 if BHP_limit_I < BHP:
829 if I.getPhase() == Well.WATER:
830 q_I_w = q_I
831 q_I_g = A_fg_I/A_fw_I * q_I
832 else:
833 q_I_g = q_I
834 q_I_w = A_fw_I/A_fg_I * q_I
835 else:
836 q_I_w = Pi_I * A_fw_I/rho_fw_I * (p_f_I - BHP_I)
837 q_I_g = Pi_I * A_fg_I/rho_fg_I * (p_f_I - BHP_I)
838 else:
839 q_I_g =0
840 q_I_w =0
841 I.setWaterProduction(q_I_w * rho_fw_I/self.rho_w.rho_0)
842 I.setGasProduction(q_I_g * rho_fg_I/self.rho_g.rho_0 )
843 return u2

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26