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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26