/[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 3485 - (show annotations)
Thu Mar 24 22:44:40 2011 UTC (8 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 25386 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.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
209 1/0
210
211
212 #=======================================================
213
214 from esys.escript.linearPDEs import LinearPDE
215 from esys.finley import Rectangle, Brick
216 from esys.escript import unitsSI as U
217 from esys.weipa import saveVTK
218 from math import pi, ceil
219 import sys
220 import time
221
222 # ======================= Default Values ==================================================
223 H=2.*U.m # height
224 L=H # length
225 NE_H=5 # number of elements in H-direction.
226 NE_L=int(ceil(L*NE_H/H))
227 DIM=2
228 T_END=40000.0*U.sec # end time at max strain = 0.5
229 DT=1.*U.sec
230
231 CHI=1.0 # =0 = explicit Euler, 0.5 = Crank-Nicholson, 1. = implicit Euler
232
233 G=-9.81*U.m/U.sec**2
234 S_W_MIN = 0.01
235 S_W_MAX = 0.75
236
237 M=3./2. *0+1
238 GAMMA=1.e-4*(1./U.Pa)
239
240 V_L=0.015*U.m**3/U.kg *0
241 P_L=6.109*U.Mega*U.Pa
242 P_A=101.325*U.Kilo*U.Pa
243 RHO_AG=0.717*U.kg/U.m**3
244 P_0G=P_A
245 RHO_0G=RHO_AG
246 RHO_W = 1000.*U.kg/U.m**3
247 RHO_COAL=1250*U.kg/U.m**3
248 ETA_W=1.e-3*U.kg/U.m/U.sec
249 ETA_G=1.1e-5*U.kg/U.m/U.sec
250 LAM=1./3.
251 PHI_0=0.00804
252 EPS_L=0.02295
253 K_S=8425*U.Mega*U.Pa
254 ALPHA=0.005*1./U.Pa *0
255 EPS_V=0
256 K_0=3.7996e-17*U.m**3
257 D_G=20e-10*0
258
259
260 TOL=1.e-4
261 OUT_DIR="ReservoirResults"
262
263
264 class Porosity(object):
265 def __init__(self, phi0=1., alpha=0., Ks=1., eps_v=0, p0=0., eps_L=0, p_L=0):
266 self.phi0=phi0
267 self.alpha=alpha
268 self.eps_v=eps_v
269 self.p_0=p0
270 self.p_L=p_L
271 self.eps_L = eps_L
272 self.Ks = Ks
273
274 def getValue(self,p=0):
275 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)))
276
277 def __call__(self, p=0):
278 return self.getValue(p)
279
280 def getDp(self, p=0):
281 return self.alpha * (1./self.Ks - self.eps_L * self.p_L/(p-self.p_L)**2 )
282
283 class Density(object):
284 def __init__(self, p_a =0, rho_a=0, p0 =0, rho0=0):
285 self.p_a=p_a
286 self.rho_a=rho_a
287 self.p0=p0
288 self.rho0=rho0
289 def getValue(self,p=0):
290 return (self.rho_a - self.rho0)/(self.p_a - self.p0)*(p-self.p0) + self.rho0
291 def __call__(self, p=0):
292 return self.getValue(p)
293
294 def getDp(self, p=0):
295 return (self.rho_a - self.rho0)/(self.p_a - self.p0)
296
297
298 class GasStorageCapacity(object):
299 def __init__(self):
300 pass
301
302 def __call__(self, p=0):
303 return self.getValue(p)
304
305 def getValue(self,p=0):
306 return 0.*p
307 def getDp(self, p=0):
308 return 0.*p
309
310 class GasStorageCapacity(object):
311 def __init__(self):
312 raise NotImplementedError
313
314 def getValue(self,p=0):
315 raise NotImplementedError
316 def getDp(self, p=0):
317 raise NotImplementedError
318 def __call__(self, p=0):
319 return self.getValue(p)
320
321
322
323 class LangmuirIsotherm(GasStorageCapacity):
324 def __init__(self, V_L=0, p_L=0., rho_a=1. ):
325 self.V_L=V_L
326 self.p_L=p_L
327 self.rho_a=rho_a
328
329 def getValue(self,p=0):
330 return (self.V_L*self.rho_a)*p/(self.p_L +p)
331
332 def getDp(self, p=0):
333 return (self.V_L*self.rho_a)*p/(self.p_L +p)**2
334
335
336 class WaterVolumeFraction(object):
337 """
338 class defining water volume fraction S_w as function of suction pressure p
339 """
340 def __init__(self):
341 raise NotImplementedError
342 def __call__(self, p=0):
343 return self.getValue(p)
344 def getValue(self,p=0):
345 """
346 returns water volume fraction for given suction pressure p
347 """
348 raise NotImplementedError
349 def getDp(self, p=0):
350 """
351 returns derivative of water volume fraction with respect to suction pressure for given suction pressure p
352 """
353 raise NotImplementedError
354
355 class VanGenuchten(WaterVolumeFraction):
356 def __init__(self, S_min=0., S_max=1., m =1, gamma=0):
357 self.S_min = S_min
358 self.S_max = S_max
359 self.m = m
360 self.gamma = gamma
361 def getEffFraction(self,s):
362 return (s-self.S_min)/(self.S_max-self.S_min)
363 def getFraction(self,s_eff):
364 return self.S_min+s_eff*(self.S_max-self.S_min)
365
366 def getValue(self,p=0):
367 s_hat = (1+(abs(p) * self.gamma)**self.m)**(1./self.m-1)
368 return self.getFraction(s_hat)
369 def getDp(self, p=0):
370 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)
371
372 if DIM==2:
373 dom=Rectangle(NE_L,NE_H,l0=L,l1=H,order=1,optimize=True)
374 else:
375 dom=Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H,order=1,optimize=True)
376
377
378
379 S_w=VanGenuchten(S_min=S_W_MIN, S_max=S_W_MAX, m =M, gamma=GAMMA)
380 C_c=LangmuirIsotherm(V_L=V_L, p_L=P_L, rho_a=RHO_AG )
381
382 Rho_g=Density(p_a = P_A, rho_a= RHO_AG, p0 = 2*P_A, rho0=RHO_0G)
383 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)
384
385 pde=LinearPDE(dom,numEquations=3, numSolutions=3)
386 pde.getSolverOptions().setVerbosityOn()
387 pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
388
389 t=0
390 n=0
391 dt=DT
392 mkDir(OUT_DIR)
393
394 z=dom.getX()[DIM-1]
395 z_top=sup(z)
396 z_bot=inf(z)
397
398
399 suction=(RHO_W-RHO_AG)*G*(z-z_top)
400 p_g=RHO_AG*G*(z-z_top)+P_A
401 m_g=( 1-S_w(suction)) * ( Rho_g(p_g) * Phi(p_g) - RHO_COAL * C_c(p_g) )
402
403 q=pde.createCoefficient("q")
404 r=pde.createCoefficient("r")
405 # fix pressures and movable mass on the surface
406 q[0]=whereZero(z-z_top)+whereZero(z-z_bot)
407 q[1]=whereZero(z-z_top)+whereZero(z-z_bot)
408 q[2]=whereZero(z-z_top)+whereZero(z-z_bot)
409 r[0]=suction
410 r[1]=m_g
411 r[2]=p_g
412 pde.setValue(q=q, r=r)
413
414 norm_p_g=Lsup(p_g)
415 norm_suction=Lsup(suction)
416
417 print " m_g %e"%(Lsup(m_g))
418 print " p_g %e"%(norm_p_g)
419 print " suction %e"%(norm_suction)
420
421 while t< T_END:
422 suction_old=suction
423 p_g_old=p_g
424 m_g_old=m_g
425 norm_dp_g = norm_p_g
426 norm_dsuction = norm_suction
427
428 print "Time step %d at %e sec"%(n+1,t+dt)
429 niter=0
430
431 while norm_dp_g > norm_p_g * TOL or norm_dsuction > norm_suction * TOL:
432 print " Iteration step %d:"%(niter,)
433
434 s_w=S_w(suction)
435 DS_wDp=S_w.getDp(suction)
436
437 rho_g=Rho_g(p_g)
438 DRho_gDp=Rho_g.getDp(p_g)
439
440 c = C_c(p_g)
441 DC_cDp=C_c.getDp(p_g)
442
443 phi=Phi(p_g)
444 DPhiDp=Phi.getDp(p_g)
445
446 print " porosity range phi:",inf(phi), sup(phi)
447 print " water volume fraction :",inf(s_w), sup(s_w)
448
449 D=pde.createCoefficient("D")
450 A=pde.createCoefficient("A")
451 X=pde.createCoefficient("X")
452 Y=pde.createCoefficient("Y")
453
454
455 # p_g equation:
456 Hp =phi*DRho_gDp + rho_g*DPhiDp - RHO_COAL*DC_cDp
457 Homega=DS_wDp*(rho_g*phi - RHO_COAL*c)
458 D[2,0]=-Homega
459 D[2,1]=-1.
460 D[2,2]= Hp
461 Y[2]=Hp * p_g_old - m_g_old - Homega*suction_old
462
463 # revise :::
464 k=K_0*(phi/PHI_0)**3
465 # s_w_eff=S_w.getEffFraction(s_w)
466 s_w_eff=s_w
467 k_hat=sqrt(s_w_eff) * (1-(1-s_w_eff**(1./LAM))**LAM)**2
468 k_w=k*k_hat
469 k_g=k*(1-k_hat)
470
471 K00=CHI*dt*k_w/ETA_W
472 K11=CHI*dt*rho_g*D_G
473 K12=CHI*dt*rho_g*k_g/ETA_G
474
475 # suction equation:
476 Y[0]=phi*DS_wDp*p_g_old+s_w*DPhiDp*suction_old
477 X[0,:] = - dt * k_w/ETA_W * ( (1-CHI) * grad(suction+p_g) - G * RHO_W * kronecker(DIM)[DIM-1] )
478 print Y[0]
479 print X[0,:]
480 D[0,0]=phi*DS_wDp
481 D[0,2]=s_w*DPhiDp
482 for i in range(DIM):
483 A[0,i,0,i]=K00
484 A[0,i,2,i]=K00
485
486
487 # p_g equation:
488 Y[1]=m_g_old
489 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]) ))
490 print X[1,:]
491 D[1,1]=1.
492 for i in range(DIM):
493 A[1,i,1,i]=K11
494 A[1,i,2,i]=K12
495
496
497 pde.setValue(A=A, D=D, X=X, Y=Y)
498 u=pde.getSolution()
499 m_g, m_g2=u[1], m_g
500 p_g, p_g2=u[2], p_g
501 suction, suction2=u[0], suction
502
503 norm_p_g=Lsup(p_g)
504 norm_suction=Lsup(suction)
505 norm_dp_g = Lsup(p_g-p_g2)
506 norm_dsuction = Lsup(suction- suction2)
507 print inf(u[0])
508 print sup(u[0])
509
510 print " m_g correction %e/%e"%(Lsup(m_g-m_g2), Lsup(m_g))
511 print " p_g correction %e/%e"%(norm_dp_g, norm_p_g)
512 print " suction correction %e/%e"%(norm_dsuction, norm_suction)
513
514 niter+=1
515 if niter >1: 1/0
516
517
518 print "Time step %d completed."%n
519 n+=1
520 t+=dt
521
522 1/0
523
524
525
526 #===============================
527 def meanValue(arg):
528 return integrate(arg,Function(arg.getDomain()))/(H*L)
529
530 #Boundary conditions:
531 # 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).
532
533 #Variables calculated and written to an output file:
534 # time
535 # differential stress (S_33-S_11)
536 # deviatoric stress (S_33 - p)
537 # Axial and transverse strain
538 # damage and damage rate
539
540
541
542 # # material parameters
543 C_V = 2.e-5/(U.Mega*U.Pa)
544 C_D = 3/U.sec
545 C_1 = 1e-10/U.sec #1e-12/U.sec
546 C_2 = 0.02 #0.03
547 XI_0 = -0.56
548 LAME_0 = 29e9*U.Pa
549 MU_0 = 19e9*U.Pa
550 ALPHA_0 = 0.
551 ALPHA_pert = 0.6
552 RHO = 2800*U.kg/U.m**3
553 G = 10*U.m/U.sec**2 *0
554 # # BC parameters
555 SIGMA_N=+50.*U.Mega*U.Pa #confining -50 , +50 gives tension -> high damage levels
556 DIM=2
557 VMAX=-1.*U.m/U.sec/800000. #*2
558 DT_MAX=50000.*U.sec
559 DT=400.*U.sec/10000.
560 #T_END=30*DT
561 if DIM==2:
562 xc=[L/2,H/2]
563 else:
564 xc=[L/2,L/2,H/2]
565 WWW=min(H,L)*0.02 # for one-element with a_0=pert+ALPHA_0 should be ~element length (=L/NE_L)
566
567 VERBOSE=True
568 DT_VIS=T_END/100 # time distane between two visulaization files
569 DN_VIS=20 # maximum counter increment between two visulaization files
570 VIS_DIR="damage-localization" # name of the directory for vis files <<<<<<<<<<<<<<<<<<<
571 DIAGNOSE_FN="diagnose_damage.csv"
572
573 ODE_TOL=0.1
574 ODE_ITER_TOL=1.e-6
575 ODE_ITER_MAX=15
576 DEPS_MAX=0.01
577 DU_ITER_TOL=1e-4
578 DU_ITER_MAX=20
579
580 diagnose=FileWriter(DIAGNOSE_FN,append=False)
581 #===================================
582 S=0.5*XI_0*((2.*MU_0+3.*LAME_0)/(3.-XI_0**2) + LAME_0)
583 GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
584
585 class DivergenceError(Exception):
586 pass
587
588 def solveODE(u0, a, b, dt):
589 """
590 solves du/dt=a*exp(b*u) u(t=0)=u0 and return approximation at t=dt
591
592 we use backwards Euler by solving u-u0 = dt * a*exp(b*u) alpha=solveODE(alpha_old, a,b, dt)
593 with newton scheme u <- u - (u-u0 - dt * a*exp(b*u)) / (1-dt*a*b*exp(b*u))
594 """
595 u=u0.copy()
596 norm_du=1.
597 norm_u=0.
598 n=0
599 while norm_du > ODE_ITER_TOL * norm_u:
600 H=-dt*a*exp(b*u)
601 du=-(u-u0+H)/(1+b*H)
602 u+=du
603 norm_du = Lsup(du)
604 norm_u = Lsup(u)
605 n+=1
606 if n>ODE_ITER_MAX:
607 print "*** ODE iteration has given up after %s steps with correction %e"%(n, norm_du)
608 raise DivergenceError,"ODE iteration failed after %s steps."%n
609 print " ODE iteration completed after %d steps with correction %e."%(n,norm_du)
610 return u
611
612 #======================
613 t=0 # time stamp
614 n=0 # time step counter
615 dt=DT # current time step size
616 t_vis=0
617 n_vis=0
618 counter_vis=0
619 mkDir(VIS_DIR)
620 #=========================
621 #
622 # ........set up domain
623 #
624
625 BBOX=boundingBox(dom)
626 DIM=dom.getDim()
627 x=dom.getX()
628 #
629 # initial values:
630 #
631 sigma=Tensor(0.,Function(dom))
632 eps_e=Tensor(0.,Function(dom))
633 #
634 # nucleation site - alpha perturbation:
635
636 alpha=ALPHA_0+ALPHA_pert*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
637 #MU_0 = MU_0 + MU_0/100*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
638
639 #saveVTK("alpha0_test-cd_0-mu0_pert.vtu",alpha=alpha) # map initial damage.vtu <<<<<<<<<<<<<<<<<<
640
641
642 # pde.getSolverOptions().setPreconditioner(pde.getSolverOptions().AMG)
643 pde.getSolverOptions().setTolerance(DU_ITER_TOL**2)
644
645
646
647 fixed_v_mask=Vector(0,Solution(dom))
648 v0=Vector(0.,ContinuousFunction(dom))
649
650 for d in range(DIM):
651 fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
652 if d == DIM-1:
653 fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)
654 v0[d]=(x[d]-BBOX[d][0])/(BBOX[d][1]-BBOX[d][0])*VMAX
655
656 pde.setValue(Y=-G*RHO*kronecker(DIM)[DIM-1], q=fixed_v_mask)
657 du=Vector(0.,Solution(dom))
658 u=Vector(0.,Solution(dom))
659 norm_du=0.
660 deps=Tensor(0,Function(dom))
661 i_eta=0
662 #
663 # ........let the show begin:
664 #
665 k3=kronecker(DIM)
666 k3Xk3=outer(k3,k3)
667 alpha_old=alpha
668 dt_old=None
669 diagnose.write("t, -e22, e11, s00-s22, alpha_max, alpha_av\n")
670
671 while t<T_END:
672
673 print "===== start time step %d ===== "%(n+1,)
674
675 eps_e_old = eps_e
676 sigma_old = sigma
677 alpha_old, alpha_oold = alpha, alpha_old
678 du_old=du
679 converged=False
680 # start the iteration for deps on a time step: deps from the last time step is used as an initial guess:
681 while not converged:
682 iter_du=0
683 if dt_old!=None:
684 du=du_old*(dt/dt_old)
685 else:
686 du=du_old
687 norm_ddu=Lsup(du)
688 norm_u = Lsup(u)
689 deps=symmetric(grad(du))
690 sigma=sigma_old
691 # start iteration :
692 try:
693 while norm_ddu > DU_ITER_TOL * norm_u or iter_du == 0 :
694 print " start iteration step %d at time step %d:"%(iter_du,n+1)
695
696 eps_e = eps_e_old + deps-(dt/2)*i_eta*deviatoric(sigma)
697
698 I1=trace(eps_e)
699 sqI2=length(eps_e)
700 xi=safeDiv(I1,sqI2)
701 # ......update damage parameter:
702 m=wherePositive(xi-XI_0)
703 a=sqI2**2*(xi-XI_0)*(m*C_D + (1-m)* C_1)
704 b=(1-m)*(1./C_2)
705
706 alpha=solveODE(alpha_old, a,b, dt)
707 if sup(alpha) > 1.:
708 print "*** damage parameter %e > 1"%(sup(alpha), )
709 raise DivergenceError,"damage parameter %e > 1"%(sup(alpha), )
710
711 alpha_dot=(alpha-alpha_old)/dt
712 i_eta = clip(2*C_V*alpha_dot,minval=0.)
713
714 gamma=alpha*GAMMA_M
715 lame=LAME_0
716 mu=MU_0*(1-alpha)
717
718 print " alpha = [ %e, %e]"%(inf(alpha),sup(alpha))
719 print " xi = [ %e, %e]"%(inf(xi),sup(xi))
720 print " gamma = [ %e, %e]"%(inf(gamma),sup(gamma))
721
722 sigma = lame * I1 * k3 + 2* mu * eps_e - gamma * ( sqI2 * k3 + xi * eps_e )
723
724 pde.setValue(A = mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3)
725 pde.setValue(X=-sigma, y=SIGMA_N*dom.getNormal(), r=dt*v0-du)
726
727
728 ddu=pde.getSolution()
729 deps+=symmetric(grad(ddu))
730 du=du+ddu
731 norm_ddu=Lsup(ddu)
732 norm_u=Lsup(u+du)
733 print " displacement change update = %e of %e"%(norm_ddu, norm_u)
734 iter_du+=1
735 if iter_du > DU_ITER_MAX:
736 print "*** displacement iteration has given up after %s steps with rel. correction %e"%(iter_du, norm_ddu/norm_u)
737 raise DivergenceError,"displacement iteration failed after %s steps."%iter_du
738 converged=True
739 print "displacement iteration converged after %d steps (rel. increment = %e)."%(iter_du, norm_ddu/norm_u)
740
741 except DivergenceError:
742 converged=False
743 dt*=0.5
744 print "*** iteration is resumed with new time step size = %e"%(dt,)
745
746 u+=du
747 n+=1
748 t+=dt
749 diagnose.write(("%e,"*6+"\n")%(t, meanValue(-symmetric(grad(u))[DIM-1,DIM-1]),
750 meanValue(symmetric(grad(u))[0,0]),
751 meanValue(sigma[0,0]-sigma[DIM-1,DIM-1]),
752 sup(alpha),
753 meanValue(alpha))) #meanValue(alpha_dot)))
754 print "time step %s (t=%e) completed."%(n,t)
755 #
756 # ............visualization:
757 #
758 # dom.setX(dom.getX()+du)
759 if t>=t_vis or n>n_vis:
760 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
761 print "visualization file %d for time step %e generated."%(counter_vis,t)
762 counter_vis+=1
763 t_vis+=DT_VIS
764 n_vis+=DN_VIS
765 #
766 # ............control time step size:
767 #
768 # ss=sup(length(deps))
769 # if ss>0:
770 # dt_new=DEPS_MAX/ss*dt
771 # print " time step size to control strain increment %s."%(dt_new,)
772 # else:
773 dt_new=dt
774 if dt_old != None:
775 dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
776 norm_alpha=Lsup(alpha)
777 fac=Lsup(dd_alpha)
778 if norm_alpha > 0: fac*=1./norm_alpha
779 if fac>0:
780 error=fac*0.5*dt**2
781 print " estimated local error for time step size %e is %e"%(dt,error)
782 dt_new=sqrt(ODE_TOL*2/fac)
783 print " new time step size to maintain error level = %e."%(dt_new,)
784 if n<10:
785 fac=10.
786 else:
787 fac=1.3
788 dt_new=min(max(dt_new,dt/fac),dt*fac) # aviod rapid changes
789 print "new time step size %e"%dt_new
790 dt, dt_old=dt_new, dt
791 print "#### Time integartion completed #########"

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26