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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3484 - (show annotations)
Thu Mar 24 05:42:07 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 25385 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 Coal Seam gasL ECLIPSE test case
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 import *
24 from esys.escript import unitsSI as U
25 from coalgas import *
26 import time
27 from esys.finley import Rectangle
28
29 L_X=168*U.km
30 L_Y=168*U.km
31 L_Z=10*U.m
32
33 N_X=21
34 N_Y=21
35
36
37 PERM_F_X = 100 * U.mDarcy
38 PERM_F_Y = 100 * U.mDarcy
39 PERM_F_Z = 1e-4 * U.mDarcy
40 PHI_F_0 = 0.01
41 P_F_0 = 69 * U.bar
42
43 # these object correspond to the ECLIPSE input files
44 PVTW={ "p_ref" : 1000 * U.bar ,
45 "B_ref" : 0.997 ,
46 "C" : 3.084E-06 /U.bar,
47 "mu_ref" : 0.68673 * U.cPoise,
48 "C_v" : 0/U.bar
49
50 }
51 GRAVITY = { "water" : 1.0,
52 "gas" : .553 }
53 ROCK = { "p_ref" : 1000 * U.bar ,
54 "C" : 3.3E-004 * 1./U.bar }
55
56 LANGMUIR = [
57 [ 0 * U.bar , 0.00000000 ],
58 [ 100 * U.bar , 0.00213886 ],
59 [ 200 * U.bar , 0.00383259 ],
60 [ 300 * U.bar , 0.00520706 ],
61 [ 400 * U.bar , 0.00634474 ],
62 [ 500 * U.bar , 0.00730199 ],
63 [ 600 * U.bar , 0.00811857 ],
64 [ 700 * U.bar , 0.00882336 ],
65 [ 800 * U.bar , 0.00943786 ],
66 [ 900 * U.bar , 0.00997836 ],
67 [ 1000 * U.bar , 0.01045748 ],
68 [ 1200 * U.bar , 0.01126912 ] ]
69
70 PVDG = [
71 [ 14.70 * U.bar ,200.3800 , 0.012025 * U.cPoise ] ,
72 [ 20.00 * U.bar ,146.0600 , 0.012030 * U.cPoise ] ,
73 [ 25.00 * U.bar ,116.1461 , 0.012034 * U.cPoise ] ,
74 [ 30.00 * U.bar ,96.3132 , 0.012038 * U.cPoise ] ,
75 [ 35.00 * U.bar ,82.2113 , 0.012043 * U.cPoise ] ,
76 [ 49.33 * U.bar ,57.7891 , 0.012055 * U.cPoise ] ,
77 [ 59.00 * U.bar ,48.0866 , 0.012064 * U.cPoise ] ,
78 [ 69.00 * U.bar ,40.9441 , 0.012073 * U.cPoise ] ,
79 [ 75.00 * U.bar ,37.5839 , 0.012078 * U.cPoise ] ,
80 [ 83.00 * U.bar ,33.8685 , 0.012085 * U.cPoise ] ,
81 [ 90.00 * U.bar ,31.1661 , 0.012092 * U.cPoise ] ,
82 [ 95.00 * U.bar ,29.4827 , 0.012097 * U.cPoise ] ,
83 [ 100.00 * U.bar ,27.9698 , 0.012101 * U.cPoise ] ,
84 [ 105.00 * U.bar ,26.6028 , 0.012106 * U.cPoise ] ,
85 [ 118.60 * U.bar ,23.4749 , 0.012119 * U.cPoise ] ,
86 [ 120.00 * U.bar ,23.1937 , 0.012120 * U.cPoise ] ,
87 [ 140.00 * U.bar ,19.7977 , 0.012140 * U.cPoise ] ,
88 [ 153.23 * U.bar ,18.0443 , 0.012153 * U.cPoise ] ,
89 [ 160.00 * U.bar ,17.2607 , 0.012159 * U.cPoise ] ,
90 [ 170.00 * U.bar ,16.2188 , 0.012169 * U.cPoise ] ,
91 [ 187.86 * U.bar ,14.6373 , 0.012188 * U.cPoise ] ,
92 [ 222.49 * U.bar ,12.3027 , 0.012224 * U.cPoise ] ,
93 [ 257.13 * U.bar ,10.6038 , 0.012262 * U.cPoise ] ,
94 [ 291.76 * U.bar ,9.3134 , 0.012301 * U.cPoise ] ,
95 [ 326.39 * U.bar ,8.3001 , 0.012341 * U.cPoise ] ,
96 [ 361.02 * U.bar ,7.4835 , 0.012383 * U.cPoise ] ,
97 [ 395.66 * U.bar ,6.8114 , 0.012425 * U.cPoise ] ,
98 [ 430.29 * U.bar ,6.2491 , 0.012470 * U.cPoise ] ,
99 [ 464.92 * U.bar ,5.7715 , 0.012515 * U.cPoise ] ,
100 [ 499.55 * U.bar ,5.3610 , 0.012562 * U.cPoise ] ,
101 [ 534.19 * U.bar ,5.0043 , 0.012610 * U.cPoise ] ,
102 [ 568.82 * U.bar ,4.6917 , 0.012659 * U.cPoise ] ,
103 [ 603.45 * U.bar ,4.4154 , 0.012710 * U.cPoise ] ,
104 [ 638.08 * U.bar ,4.1695 , 0.012762 * U.cPoise ] ,
105 [ 672.72 * U.bar ,3.9491 , 0.012815 * U.cPoise ] ,
106 [ 707.35 * U.bar ,3.7507 , 0.012869 * U.cPoise ] ,
107 [ 741.98 * U.bar ,3.5711 , 0.012925 * U.cPoise ] ,
108 [ 776.61 * U.bar ,3.4076 , 0.012982 * U.cPoise ] ,
109 [ 811.25 * U.bar ,3.2583 , 0.013041 * U.cPoise ] ,
110 [ 845.88 * U.bar ,3.1214 , 0.013100 * U.cPoise ] ,
111 [ 880.51 * U.bar ,2.9953 , 0.013161 * U.cPoise ] ,
112 [ 915.14 * U.bar ,2.8790 , 0.013223 * U.cPoise ] ,
113 [ 949.78 * U.bar ,2.7712 , 0.013287 * U.cPoise ] ,
114 [ 984.41 * U.bar ,2.6711 , 0.013352 * U.cPoise ] ,
115 [ 1019.00 * U.bar ,2.5781 , 0.013418 * U.cPoise ] ,
116 [ 1053.70 * U.bar ,2.4909 , 0.013486 * U.cPoise ] ,
117 [ 1088.30 * U.bar ,2.4096 , 0.013554 * U.cPoise ] ,
118 [ 1122.90 * U.bar ,2.3334 , 0.013624 * U.cPoise ] ,
119 [ 1157.60 * U.bar ,2.2616 , 0.013696 * U.cPoise ] ,
120 [ 1192.20 * U.bar ,2.1942 , 0.013768 * U.cPoise ] ,
121 [ 1226.80 * U.bar ,2.1307 , 0.013842 * U.cPoise ] ,
122 [ 1261.50 * U.bar ,2.0705 , 0.013917 * U.cPoise ] ,
123 [ 1296.10 * U.bar ,2.0138 , 0.013994 * U.cPoise ] ,
124 [ 1330.70 * U.bar ,1.9600 , 0.014072 * U.cPoise ] ,
125 [ 1365.40 * U.bar ,1.9089 , 0.014151 * U.cPoise ] ]
126
127
128 SGFN = [
129 [ 0 , 0 , 0 ],
130 [ 0.05 , 0 , 0 ],
131 [ 0.1333 , 0.00610 , 0 ],
132 [ 0.2167 , 0.02990 , 0 ],
133 [ 0.3 , 0.0759 , 0 ],
134 [ 0.3833 , 0.1471 , 0 ],
135 [ 0.46667 , 0.2458 , 0 ],
136 [ 0.55 , 0.3739 , 0 ],
137 [ 0.6333 , 0.53300 , 0 ],
138 [ 0.7167 , 0.7246 , 0 ],
139 [ 0.8 , 0.95 , 0 ] ]
140 SWFN = [
141 [ 0.20000 , 0.00000, 0 ],
142 [ 0.28330 , 0.03280, 0 ],
143 [ 0.36670 , 0.09270, 0 ],
144 [ 0.45000 , 0.17030, 0 ],
145 [ 0.53330 , 0.26220, 0 ],
146 [ 0.61670 , 0.36650, 0 ],
147 [ 0.70000 , 0.48170, 0 ],
148 [ 0.78330 , 0.60710, 0 ],
149 [ 0.86670 , 0.74170, 0 ],
150 [ 0.95000 , 0.88500, 0 ],
151 [ 1.00000 , 1.00000, 0 ] ]
152
153
154 wellspecs = {
155 'P1' : { "X0" : [106, 106],
156 "r" : 0.253*U.m,
157 "s" : 0,
158 "Q" : 2000*U.Barrel/U.day*GRAVITY["water"],
159 "BHP" : 75*U.psi,
160 "schedule" : [0.*U.yr, 4*U.yr]
161 }
162 }
163
164 CELL_X=L_X/N_X
165 CELL_Y=L_Y/N_Y
166 # print input
167 print("<%s> Execution started."%time.asctime())
168
169 print "length x-direction = %f km"%(L_X/U.km)
170 print "number of cells in x direction = %d"%N_X
171 print "cell size in x direction = %f m"%(CELL_X/U.m)
172 print "length y-direction = %f km"%(L_Y/U.km)
173 print "number of cells in y direction = %d"%N_Y
174 print "cell size in y direction = %f m"%(CELL_Y/U.m)
175 print "fracture permeability in x direction= %f mD"%(PERM_F_X/(U.mDarcy))
176 print "fracture permeability in y direction= %f mD"%(PERM_F_Y/(U.mDarcy))
177 print "fracture permeability in z direction= %f mD"%(PERM_F_Z/(U.mDarcy))
178 print "initial porosity in fractured rock= %f"%PHI_F_0
179
180 domain=Rectangle(N_X, N_Y, l0=L_X, l1=L_Y)
181 print("<%s> Mesh set up completed."%time.asctime())
182
183
184 model = PorosityOneHalfModel(domain,
185 phi_f=Porosity(phi_0=PHI_F_0, p_0=P_F_0, p_ref=ROCK["p_ref"], C = ROCK["C"]),
186 L_g=InterpolationTable([ l[0] for l in LANGMUIR ], [ l[1] for l in LANGMUIR ] ),
187 perm_f_0=PERM_F_X,
188 perm_f_1=PERM_F_Y,
189 perm_f_2=PERM_F_Z,
190 k_w =InterpolationTable([ l[0] for l in SWFN ], [ l[1] for l in SWFN ] ),
191 k_g= InterpolationTable([ l[0] for l in SGFN ], [ l[1] for l in SGFN ] ),
192 mu_w = WaterViscosity(mu_ref = PVTW["mu_ref"], p_ref=PVTW["p_ref"], C=PVTW["C_v"]),
193 mu_g = InterpolationTable([ l[0] for l in PVDG ], [ l[2] for l in PVDG ] ),
194 rho_w = WaterDensity(B_ref=PVTW["B_ref"], p_ref = PVTW["p_ref"], C=PVTW["C"], gravity=GRAVITY["water"]),
195 rho_g=GasDensity( p = [ l[0] for l in PVDG ], B = [ l[1] for l in PVDG ], gravity=GRAVITY["gas"]),
196 wells=[ VerticalPeacemanWell(i,BHP=wellspecs[i]["BHP"],
197 Q=wellspecs[i]["Q"],
198 r=wellspecs[i]["r"],
199 X0=[ (wellspecs[i]["X0"][0]+0.5)*CELL_X, (wellspecs[i]["X0"][1]+0.5)*CELL_Y],
200 D=[CELL_X, CELL_Y, L_Z],
201 perm=[PERM_F_X, PERM_F_Y, PERM_F_Z],
202 schedule=wellspecs[i]["schedule"],
203 s=wellspecs[i]["s"]) for i in wellspecs]
204 )
205
206 model.setInitialState(p=P_F_0, S_fg=0, C_mg=None)
207
208 1/0
209
210
211 #=======================================================
212
213 from esys.escript.linearPDEs import LinearPDE
214 from esys.finley import Rectangle, Brick
215 from esys.escript import unitsSI as U
216 from esys.weipa import saveVTK
217 from math import pi, ceil
218 import sys
219 import time
220
221 # ======================= Default Values ==================================================
222 H=2.*U.m # height
223 L=H # length
224 NE_H=5 # number of elements in H-direction.
225 NE_L=int(ceil(L*NE_H/H))
226 DIM=2
227 T_END=40000.0*U.sec # end time at max strain = 0.5
228 DT=1.*U.sec
229
230 CHI=1.0 # =0 = explicit Euler, 0.5 = Crank-Nicholson, 1. = implicit Euler
231
232 G=-9.81*U.m/U.sec**2
233 S_W_MIN = 0.01
234 S_W_MAX = 0.75
235
236 M=3./2. *0+1
237 GAMMA=1.e-4*(1./U.Pa)
238
239 V_L=0.015*U.m**3/U.kg *0
240 P_L=6.109*U.Mega*U.Pa
241 P_A=101.325*U.Kilo*U.Pa
242 RHO_AG=0.717*U.kg/U.m**3
243 P_0G=P_A
244 RHO_0G=RHO_AG
245 RHO_W = 1000.*U.kg/U.m**3
246 RHO_COAL=1250*U.kg/U.m**3
247 ETA_W=1.e-3*U.kg/U.m/U.sec
248 ETA_G=1.1e-5*U.kg/U.m/U.sec
249 LAM=1./3.
250 PHI_0=0.00804
251 EPS_L=0.02295
252 K_S=8425*U.Mega*U.Pa
253 ALPHA=0.005*1./U.Pa *0
254 EPS_V=0
255 K_0=3.7996e-17*U.m**3
256 D_G=20e-10*0
257
258
259 TOL=1.e-4
260 OUT_DIR="ReservoirResults"
261
262
263 class Porosity(object):
264 def __init__(self, phi0=1., alpha=0., Ks=1., eps_v=0, p0=0., eps_L=0, p_L=0):
265 self.phi0=phi0
266 self.alpha=alpha
267 self.eps_v=eps_v
268 self.p_0=p0
269 self.p_L=p_L
270 self.eps_L = eps_L
271 self.Ks = Ks
272
273 def getValue(self,p=0):
274 return self.phi0+ self.alpha*(self.eps_v + (p-self.p_0)/self.Ks + self.eps_L * self.p_L * (self.p_0-p)/((self.p_L + self.p_0) * (p + self.p_L)))
275
276 def __call__(self, p=0):
277 return self.getValue(p)
278
279 def getDp(self, p=0):
280 return self.alpha * (1./self.Ks - self.eps_L * self.p_L/(p-self.p_L)**2 )
281
282 class Density(object):
283 def __init__(self, p_a =0, rho_a=0, p0 =0, rho0=0):
284 self.p_a=p_a
285 self.rho_a=rho_a
286 self.p0=p0
287 self.rho0=rho0
288 def getValue(self,p=0):
289 return (self.rho_a - self.rho0)/(self.p_a - self.p0)*(p-self.p0) + self.rho0
290 def __call__(self, p=0):
291 return self.getValue(p)
292
293 def getDp(self, p=0):
294 return (self.rho_a - self.rho0)/(self.p_a - self.p0)
295
296
297 class GasStorageCapacity(object):
298 def __init__(self):
299 pass
300
301 def __call__(self, p=0):
302 return self.getValue(p)
303
304 def getValue(self,p=0):
305 return 0.*p
306 def getDp(self, p=0):
307 return 0.*p
308
309 class GasStorageCapacity(object):
310 def __init__(self):
311 raise NotImplementedError
312
313 def getValue(self,p=0):
314 raise NotImplementedError
315 def getDp(self, p=0):
316 raise NotImplementedError
317 def __call__(self, p=0):
318 return self.getValue(p)
319
320
321
322 class LangmuirIsotherm(GasStorageCapacity):
323 def __init__(self, V_L=0, p_L=0., rho_a=1. ):
324 self.V_L=V_L
325 self.p_L=p_L
326 self.rho_a=rho_a
327
328 def getValue(self,p=0):
329 return (self.V_L*self.rho_a)*p/(self.p_L +p)
330
331 def getDp(self, p=0):
332 return (self.V_L*self.rho_a)*p/(self.p_L +p)**2
333
334
335 class WaterVolumeFraction(object):
336 """
337 class defining water volume fraction S_w as function of suction pressure p
338 """
339 def __init__(self):
340 raise NotImplementedError
341 def __call__(self, p=0):
342 return self.getValue(p)
343 def getValue(self,p=0):
344 """
345 returns water volume fraction for given suction pressure p
346 """
347 raise NotImplementedError
348 def getDp(self, p=0):
349 """
350 returns derivative of water volume fraction with respect to suction pressure for given suction pressure p
351 """
352 raise NotImplementedError
353
354 class VanGenuchten(WaterVolumeFraction):
355 def __init__(self, S_min=0., S_max=1., m =1, gamma=0):
356 self.S_min = S_min
357 self.S_max = S_max
358 self.m = m
359 self.gamma = gamma
360 def getEffFraction(self,s):
361 return (s-self.S_min)/(self.S_max-self.S_min)
362 def getFraction(self,s_eff):
363 return self.S_min+s_eff*(self.S_max-self.S_min)
364
365 def getValue(self,p=0):
366 s_hat = (1+(abs(p) * self.gamma)**self.m)**(1./self.m-1)
367 return self.getFraction(s_hat)
368 def getDp(self, p=0):
369 return (self.S_max-self.S_min) * (self.m-1)/self.m*(1+(self.gamma * abs(p))**self.m)**(1./self.m) * self.gamma**self.m * abs(p)**(self.m-1)*sign(p)
370
371 if DIM==2:
372 dom=Rectangle(NE_L,NE_H,l0=L,l1=H,order=1,optimize=True)
373 else:
374 dom=Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H,order=1,optimize=True)
375
376
377
378 S_w=VanGenuchten(S_min=S_W_MIN, S_max=S_W_MAX, m =M, gamma=GAMMA)
379 C_c=LangmuirIsotherm(V_L=V_L, p_L=P_L, rho_a=RHO_AG )
380
381 Rho_g=Density(p_a = P_A, rho_a= RHO_AG, p0 = 2*P_A, rho0=RHO_0G)
382 Phi=Porosity(phi0=PHI_0, alpha=ALPHA, Ks=K_S, eps_v=EPS_V, p0=P_A, eps_L=EPS_L, p_L=P_L)
383
384 pde=LinearPDE(dom,numEquations=3, numSolutions=3)
385 pde.getSolverOptions().setVerbosityOn()
386 pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
387
388 t=0
389 n=0
390 dt=DT
391 mkDir(OUT_DIR)
392
393 z=dom.getX()[DIM-1]
394 z_top=sup(z)
395 z_bot=inf(z)
396
397
398 suction=(RHO_W-RHO_AG)*G*(z-z_top)
399 p_g=RHO_AG*G*(z-z_top)+P_A
400 m_g=( 1-S_w(suction)) * ( Rho_g(p_g) * Phi(p_g) - RHO_COAL * C_c(p_g) )
401
402 q=pde.createCoefficient("q")
403 r=pde.createCoefficient("r")
404 # fix pressures and movable mass on the surface
405 q[0]=whereZero(z-z_top)+whereZero(z-z_bot)
406 q[1]=whereZero(z-z_top)+whereZero(z-z_bot)
407 q[2]=whereZero(z-z_top)+whereZero(z-z_bot)
408 r[0]=suction
409 r[1]=m_g
410 r[2]=p_g
411 pde.setValue(q=q, r=r)
412
413 norm_p_g=Lsup(p_g)
414 norm_suction=Lsup(suction)
415
416 print " m_g %e"%(Lsup(m_g))
417 print " p_g %e"%(norm_p_g)
418 print " suction %e"%(norm_suction)
419
420 while t< T_END:
421 suction_old=suction
422 p_g_old=p_g
423 m_g_old=m_g
424 norm_dp_g = norm_p_g
425 norm_dsuction = norm_suction
426
427 print "Time step %d at %e sec"%(n+1,t+dt)
428 niter=0
429
430 while norm_dp_g > norm_p_g * TOL or norm_dsuction > norm_suction * TOL:
431 print " Iteration step %d:"%(niter,)
432
433 s_w=S_w(suction)
434 DS_wDp=S_w.getDp(suction)
435
436 rho_g=Rho_g(p_g)
437 DRho_gDp=Rho_g.getDp(p_g)
438
439 c = C_c(p_g)
440 DC_cDp=C_c.getDp(p_g)
441
442 phi=Phi(p_g)
443 DPhiDp=Phi.getDp(p_g)
444
445 print " porosity range phi:",inf(phi), sup(phi)
446 print " water volume fraction :",inf(s_w), sup(s_w)
447
448 D=pde.createCoefficient("D")
449 A=pde.createCoefficient("A")
450 X=pde.createCoefficient("X")
451 Y=pde.createCoefficient("Y")
452
453
454 # p_g equation:
455 Hp =phi*DRho_gDp + rho_g*DPhiDp - RHO_COAL*DC_cDp
456 Homega=DS_wDp*(rho_g*phi - RHO_COAL*c)
457 D[2,0]=-Homega
458 D[2,1]=-1.
459 D[2,2]= Hp
460 Y[2]=Hp * p_g_old - m_g_old - Homega*suction_old
461
462 # revise :::
463 k=K_0*(phi/PHI_0)**3
464 # s_w_eff=S_w.getEffFraction(s_w)
465 s_w_eff=s_w
466 k_hat=sqrt(s_w_eff) * (1-(1-s_w_eff**(1./LAM))**LAM)**2
467 k_w=k*k_hat
468 k_g=k*(1-k_hat)
469
470 K00=CHI*dt*k_w/ETA_W
471 K11=CHI*dt*rho_g*D_G
472 K12=CHI*dt*rho_g*k_g/ETA_G
473
474 # suction equation:
475 Y[0]=phi*DS_wDp*p_g_old+s_w*DPhiDp*suction_old
476 X[0,:] = - dt * k_w/ETA_W * ( (1-CHI) * grad(suction+p_g) - G * RHO_W * kronecker(DIM)[DIM-1] )
477 print Y[0]
478 print X[0,:]
479 D[0,0]=phi*DS_wDp
480 D[0,2]=s_w*DPhiDp
481 for i in range(DIM):
482 A[0,i,0,i]=K00
483 A[0,i,2,i]=K00
484
485
486 # p_g equation:
487 Y[1]=m_g_old
488 X[1,:] = - dt * rho_g * ((1-CHI) * ( D_G * grad(m_g_old) + rho_g * k_g/ETA_G * ( grad(p_g_old) - G * rho_g * kronecker(DIM)[DIM-1]) ))
489 print X[1,:]
490 D[1,1]=1.
491 for i in range(DIM):
492 A[1,i,1,i]=K11
493 A[1,i,2,i]=K12
494
495
496 pde.setValue(A=A, D=D, X=X, Y=Y)
497 u=pde.getSolution()
498 m_g, m_g2=u[1], m_g
499 p_g, p_g2=u[2], p_g
500 suction, suction2=u[0], suction
501
502 norm_p_g=Lsup(p_g)
503 norm_suction=Lsup(suction)
504 norm_dp_g = Lsup(p_g-p_g2)
505 norm_dsuction = Lsup(suction- suction2)
506 print inf(u[0])
507 print sup(u[0])
508
509 print " m_g correction %e/%e"%(Lsup(m_g-m_g2), Lsup(m_g))
510 print " p_g correction %e/%e"%(norm_dp_g, norm_p_g)
511 print " suction correction %e/%e"%(norm_dsuction, norm_suction)
512
513 niter+=1
514 if niter >1: 1/0
515
516
517 print "Time step %d completed."%n
518 n+=1
519 t+=dt
520
521 1/0
522
523
524
525 #===============================
526 def meanValue(arg):
527 return integrate(arg,Function(arg.getDomain()))/(H*L)
528
529 #Boundary conditions:
530 # axial loading: they applied a stress inversely proportional to the acoustic emission rate. We could have the axial forcing a stress or velocity inversely proportional to dalpha/dt (only when it is positive, and with the applied forcing rate going to zero when damage accumulation rate goes to a value we can determine in a test run with constant forcing). If this is to challenging or time consuming we could have a constant axial strain rate with very short time steps (at least when alpha increases above 0.3).
531
532 #Variables calculated and written to an output file:
533 # time
534 # differential stress (S_33-S_11)
535 # deviatoric stress (S_33 - p)
536 # Axial and transverse strain
537 # damage and damage rate
538
539
540
541 # # material parameters
542 C_V = 2.e-5/(U.Mega*U.Pa)
543 C_D = 3/U.sec
544 C_1 = 1e-10/U.sec #1e-12/U.sec
545 C_2 = 0.02 #0.03
546 XI_0 = -0.56
547 LAME_0 = 29e9*U.Pa
548 MU_0 = 19e9*U.Pa
549 ALPHA_0 = 0.
550 ALPHA_pert = 0.6
551 RHO = 2800*U.kg/U.m**3
552 G = 10*U.m/U.sec**2 *0
553 # # BC parameters
554 SIGMA_N=+50.*U.Mega*U.Pa #confining -50 , +50 gives tension -> high damage levels
555 DIM=2
556 VMAX=-1.*U.m/U.sec/800000. #*2
557 DT_MAX=50000.*U.sec
558 DT=400.*U.sec/10000.
559 #T_END=30*DT
560 if DIM==2:
561 xc=[L/2,H/2]
562 else:
563 xc=[L/2,L/2,H/2]
564 WWW=min(H,L)*0.02 # for one-element with a_0=pert+ALPHA_0 should be ~element length (=L/NE_L)
565
566 VERBOSE=True
567 DT_VIS=T_END/100 # time distane between two visulaization files
568 DN_VIS=20 # maximum counter increment between two visulaization files
569 VIS_DIR="damage-localization" # name of the directory for vis files <<<<<<<<<<<<<<<<<<<
570 DIAGNOSE_FN="diagnose_damage.csv"
571
572 ODE_TOL=0.1
573 ODE_ITER_TOL=1.e-6
574 ODE_ITER_MAX=15
575 DEPS_MAX=0.01
576 DU_ITER_TOL=1e-4
577 DU_ITER_MAX=20
578
579 diagnose=FileWriter(DIAGNOSE_FN,append=False)
580 #===================================
581 S=0.5*XI_0*((2.*MU_0+3.*LAME_0)/(3.-XI_0**2) + LAME_0)
582 GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
583
584 class DivergenceError(Exception):
585 pass
586
587 def solveODE(u0, a, b, dt):
588 """
589 solves du/dt=a*exp(b*u) u(t=0)=u0 and return approximation at t=dt
590
591 we use backwards Euler by solving u-u0 = dt * a*exp(b*u) alpha=solveODE(alpha_old, a,b, dt)
592 with newton scheme u <- u - (u-u0 - dt * a*exp(b*u)) / (1-dt*a*b*exp(b*u))
593 """
594 u=u0.copy()
595 norm_du=1.
596 norm_u=0.
597 n=0
598 while norm_du > ODE_ITER_TOL * norm_u:
599 H=-dt*a*exp(b*u)
600 du=-(u-u0+H)/(1+b*H)
601 u+=du
602 norm_du = Lsup(du)
603 norm_u = Lsup(u)
604 n+=1
605 if n>ODE_ITER_MAX:
606 print "*** ODE iteration has given up after %s steps with correction %e"%(n, norm_du)
607 raise DivergenceError,"ODE iteration failed after %s steps."%n
608 print " ODE iteration completed after %d steps with correction %e."%(n,norm_du)
609 return u
610
611 #======================
612 t=0 # time stamp
613 n=0 # time step counter
614 dt=DT # current time step size
615 t_vis=0
616 n_vis=0
617 counter_vis=0
618 mkDir(VIS_DIR)
619 #=========================
620 #
621 # ........set up domain
622 #
623
624 BBOX=boundingBox(dom)
625 DIM=dom.getDim()
626 x=dom.getX()
627 #
628 # initial values:
629 #
630 sigma=Tensor(0.,Function(dom))
631 eps_e=Tensor(0.,Function(dom))
632 #
633 # nucleation site - alpha perturbation:
634
635 alpha=ALPHA_0+ALPHA_pert*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
636 #MU_0 = MU_0 + MU_0/100*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
637
638 #saveVTK("alpha0_test-cd_0-mu0_pert.vtu",alpha=alpha) # map initial damage.vtu <<<<<<<<<<<<<<<<<<
639
640
641 # pde.getSolverOptions().setPreconditioner(pde.getSolverOptions().AMG)
642 pde.getSolverOptions().setTolerance(DU_ITER_TOL**2)
643
644
645
646 fixed_v_mask=Vector(0,Solution(dom))
647 v0=Vector(0.,ContinuousFunction(dom))
648
649 for d in range(DIM):
650 fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
651 if d == DIM-1:
652 fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)
653 v0[d]=(x[d]-BBOX[d][0])/(BBOX[d][1]-BBOX[d][0])*VMAX
654
655 pde.setValue(Y=-G*RHO*kronecker(DIM)[DIM-1], q=fixed_v_mask)
656 du=Vector(0.,Solution(dom))
657 u=Vector(0.,Solution(dom))
658 norm_du=0.
659 deps=Tensor(0,Function(dom))
660 i_eta=0
661 #
662 # ........let the show begin:
663 #
664 k3=kronecker(DIM)
665 k3Xk3=outer(k3,k3)
666 alpha_old=alpha
667 dt_old=None
668 diagnose.write("t, -e22, e11, s00-s22, alpha_max, alpha_av\n")
669
670 while t<T_END:
671
672 print "===== start time step %d ===== "%(n+1,)
673
674 eps_e_old = eps_e
675 sigma_old = sigma
676 alpha_old, alpha_oold = alpha, alpha_old
677 du_old=du
678 converged=False
679 # start the iteration for deps on a time step: deps from the last time step is used as an initial guess:
680 while not converged:
681 iter_du=0
682 if dt_old!=None:
683 du=du_old*(dt/dt_old)
684 else:
685 du=du_old
686 norm_ddu=Lsup(du)
687 norm_u = Lsup(u)
688 deps=symmetric(grad(du))
689 sigma=sigma_old
690 # start iteration :
691 try:
692 while norm_ddu > DU_ITER_TOL * norm_u or iter_du == 0 :
693 print " start iteration step %d at time step %d:"%(iter_du,n+1)
694
695 eps_e = eps_e_old + deps-(dt/2)*i_eta*deviatoric(sigma)
696
697 I1=trace(eps_e)
698 sqI2=length(eps_e)
699 xi=safeDiv(I1,sqI2)
700 # ......update damage parameter:
701 m=wherePositive(xi-XI_0)
702 a=sqI2**2*(xi-XI_0)*(m*C_D + (1-m)* C_1)
703 b=(1-m)*(1./C_2)
704
705 alpha=solveODE(alpha_old, a,b, dt)
706 if sup(alpha) > 1.:
707 print "*** damage parameter %e > 1"%(sup(alpha), )
708 raise DivergenceError,"damage parameter %e > 1"%(sup(alpha), )
709
710 alpha_dot=(alpha-alpha_old)/dt
711 i_eta = clip(2*C_V*alpha_dot,minval=0.)
712
713 gamma=alpha*GAMMA_M
714 lame=LAME_0
715 mu=MU_0*(1-alpha)
716
717 print " alpha = [ %e, %e]"%(inf(alpha),sup(alpha))
718 print " xi = [ %e, %e]"%(inf(xi),sup(xi))
719 print " gamma = [ %e, %e]"%(inf(gamma),sup(gamma))
720
721 sigma = lame * I1 * k3 + 2* mu * eps_e - gamma * ( sqI2 * k3 + xi * eps_e )
722
723 pde.setValue(A = mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3)
724 pde.setValue(X=-sigma, y=SIGMA_N*dom.getNormal(), r=dt*v0-du)
725
726
727 ddu=pde.getSolution()
728 deps+=symmetric(grad(ddu))
729 du=du+ddu
730 norm_ddu=Lsup(ddu)
731 norm_u=Lsup(u+du)
732 print " displacement change update = %e of %e"%(norm_ddu, norm_u)
733 iter_du+=1
734 if iter_du > DU_ITER_MAX:
735 print "*** displacement iteration has given up after %s steps with rel. correction %e"%(iter_du, norm_ddu/norm_u)
736 raise DivergenceError,"displacement iteration failed after %s steps."%iter_du
737 converged=True
738 print "displacement iteration converged after %d steps (rel. increment = %e)."%(iter_du, norm_ddu/norm_u)
739
740 except DivergenceError:
741 converged=False
742 dt*=0.5
743 print "*** iteration is resumed with new time step size = %e"%(dt,)
744
745 u+=du
746 n+=1
747 t+=dt
748 diagnose.write(("%e,"*6+"\n")%(t, meanValue(-symmetric(grad(u))[DIM-1,DIM-1]),
749 meanValue(symmetric(grad(u))[0,0]),
750 meanValue(sigma[0,0]-sigma[DIM-1,DIM-1]),
751 sup(alpha),
752 meanValue(alpha))) #meanValue(alpha_dot)))
753 print "time step %s (t=%e) completed."%(n,t)
754 #
755 # ............visualization:
756 #
757 # dom.setX(dom.getX()+du)
758 if t>=t_vis or n>n_vis:
759 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),u=u, alpha=alpha, I1=trace(eps_e), eps_12=(symmetric(grad(u))[0,1]), xi=safeDiv(trace(eps_e),length(eps_e)), alpha_dot=alpha_dot, sqI2=sqI2) # tau=sigma[0,1], I2=length(eps_e)**2
760 print "visualization file %d for time step %e generated."%(counter_vis,t)
761 counter_vis+=1
762 t_vis+=DT_VIS
763 n_vis+=DN_VIS
764 #
765 # ............control time step size:
766 #
767 # ss=sup(length(deps))
768 # if ss>0:
769 # dt_new=DEPS_MAX/ss*dt
770 # print " time step size to control strain increment %s."%(dt_new,)
771 # else:
772 dt_new=dt
773 if dt_old != None:
774 dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
775 norm_alpha=Lsup(alpha)
776 fac=Lsup(dd_alpha)
777 if norm_alpha > 0: fac*=1./norm_alpha
778 if fac>0:
779 error=fac*0.5*dt**2
780 print " estimated local error for time step size %e is %e"%(dt,error)
781 dt_new=sqrt(ODE_TOL*2/fac)
782 print " new time step size to maintain error level = %e."%(dt_new,)
783 if n<10:
784 fac=10.
785 else:
786 fac=1.3
787 dt_new=min(max(dt_new,dt/fac),dt*fac) # aviod rapid changes
788 print "new time step size %e"%dt_new
789 dt, dt_old=dt_new, dt
790 print "#### Time integartion completed #########"

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26