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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3289 - (show annotations)
Wed Oct 20 07:07:11 2010 UTC (8 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 10769 byte(s)
missing factor fixed.
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 Damage mechanics
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.linearPDEs import LinearPDESystem
25 from esys.finley import Rectangle, Brick
26 from esys.escript import unitsSI as U
27 from math import pi, ceil
28 import sys
29 import time
30
31 # ======================= Default Values ==================================================
32 H=1.*U.m # height
33 L=1.*H # length
34 NE_H=1 # number of elements in H-direction.
35 NE_L=int(ceil(L*NE_H/H))
36 CASE=1
37
38 #Boundary conditions:
39 # 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).
40
41 #Variables calculated and written to an output file:
42 # time
43 # differential stress (S_33-S_11)
44 # deviatoric stress (S_33 - p)
45 # Axial and transverse strain
46 # damage and damage rate
47
48
49 T_END=60000000.0*U.sec # end time
50
51 if CASE==1 or CASE==2:
52 C_V = 2.e-5/(U.Mega*U.Pa)
53 C_D = 3/U.sec
54 C_1 = 1e-12/U.sec
55 C_2 = 0.03
56 XI_0 = -0.56
57 LAME_0 = 29e9*U.Pa
58 MU_0 = 19e9*U.Pa
59 if CASE==2:
60 ALPHA_0 = 0.999
61 else:
62 ALPHA_0 = 0.0
63 RHO = 2800*U.kg/U.m**3
64 G = 10*U.m/U.sec**2 *0
65 SIGMA_N=50.*U.Mega*U.Pa
66 DIM=3
67 VMAX=-1.*U.m/U.sec/166667.
68 DT_MAX=50.*U.sec
69 DT=DT_MAX/1000000.
70 xc=[L/2,L/2,H/2]
71 WWW=min(H,L)*0.01
72
73 else:
74 C_V = 3e-11/U.Pa
75 C_D = 5/U.sec
76 C_1 = 1e-12/U.sec
77 C_2 = 0.03
78 XI_0 = -0.8
79 LAME_0 = 46e9*U.Pa
80 MU_0 = 30e9*U.Pa
81 if CASE==3:
82 ALPHA_0 = 0.01
83 else:
84 ALPHA_0 = 0.0
85 RHO = 2800*U.kg/U.m**3
86 G = 10*U.m/U.sec**2 *0
87 VMAX=-1*U.m/U.sec
88 DT_MAX=500.*U.sec
89 DT=DT_MAX/100000.
90 SIGMA_N=0
91 DIM=2
92 xc=[L/2,H/2]
93 WWW=min(H,L)*0.08
94
95 VERBOSE=True
96 DT_VIS=T_END/100 # time distane between two visulaization files
97 DN_VIS=1 # maximum counter increment between two visulaization files
98 VIS_DIR="results" # name of the director for vis files
99
100 ODE_TOL=0.01
101 ODE_ITER_TOL=1.e-8
102 ODE_ITER_MAX=15
103 DEPS_MAX=0.01
104 TOL_DU=1e-5
105 UPDATE_OPERATOR = False
106
107
108 diagnose=FileWriter("diagnose1.csv",append=False)
109 #===================================
110 S=0.5*XI_0*((2.*MU_0+3.*LAME_0)/(3.-XI_0**2) + LAME_0)
111 GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
112 def solveODE(u0, a, b, dt):
113 """
114 solves du/dt=a*exp(b*u) u(t=0)=u0 and return approximation at t=dt
115
116 we sue backwards Euler by solving u-u0 = dt * a*exp(b*u)
117 with newton scheme u <- u - (u-u0 - dt * a*exp(b*u)) / (1-dt*a*b*exp(b*u))
118 """
119 u=u0.copy()
120 norm_du=1.
121 norm_u=0.
122 n=0
123 while norm_du > ODE_ITER_TOL * norm_u:
124 H=-dt*a*exp(b*u)
125 du=-(u-u0+H)/(1+b*H)
126 u+=du
127 norm_du = Lsup(du)
128 norm_u = Lsup(u)
129 n+=1
130 if n>ODE_ITER_MAX: raise ValueError,"ODE iteration failed."
131 print "\tODE iteration completed after %d steps."%(n,)
132 return u
133
134 #======================
135 t=0 # time stamp
136 n=0 # time step counter
137 dt=DT # current time step size
138 t_vis=0
139 n_vis=0
140 counter_vis=0
141 mkDir(VIS_DIR)
142 #=========================
143 #
144 # set up domain
145 #
146 if DIM==2:
147 dom=Rectangle(NE_L,NE_H,l0=L,l1=H,order=-1,optimize=True)
148 else:
149 dom=Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H,order=-1,optimize=True)
150
151 BBOX=boundingBox(dom)
152 DIM=dom.getDim()
153 x=dom.getX()
154 #
155 # initial values:
156 #
157 sigma=Tensor(0.,Function(dom))
158 eps_e=Tensor(0.,Function(dom))
159 if CASE==2 or CASE==3:
160 alpha=ALPHA_0*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
161 else:
162 alpha=Scalar(ALPHA_0,Function(dom))
163
164 pde=LinearPDESystem(dom)
165 pde.setSymmetryOn()
166 pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
167
168
169
170 fixed_v_mask=Vector(0,Solution(dom))
171 v0=Vector(0.,ContinuousFunction(dom))
172 if CASE == 1 or CASE==2:
173 for d in range(DIM):
174 fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
175 if d == DIM-1:
176 fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)
177 v0[d]=(x[d]-BBOX[d][0])/(BBOX[d][1]-BBOX[d][0])*VMAX
178 else:
179 for d in range(DIM):
180 fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
181 if d == 0:
182 fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)
183 v0[d]=(x[d]-BBOX[d][0])/(BBOX[d][1]-BBOX[d][0])*VMAX
184
185 pde.setValue(Y=-G*RHO*kronecker(DIM)[DIM-1], q=fixed_v_mask)
186 du=Vector(0.,Solution(dom))
187 norm_du=0.
188 deps=Tensor(0,Function(dom))
189 i_eta=0
190 #
191 # let the show begin:
192 #
193 k3=kronecker(DIM)
194 k3Xk3=outer(k3,k3)
195 alpha_old=alpha
196 dt_old=None
197 diagnose.write("t, -e22, e11, s00-s22, mu_eff, lame_eff, xi, gamma, alpha, alpha_dot\n")
198
199 while t<T_END:
200
201 print "start time step %d"%(n+1,)
202
203 eps_e_old = eps_e
204 sigma_old = sigma
205 alpha_old, alpha_oold = alpha, alpha_old
206 # start the iteration for deps on a time step: deps from the last time step is used as an initial guess:
207 iter=0
208 norm_ddu=norm_du
209 while norm_ddu > TOL_DU * norm_du or iter == 0 :
210
211 print "\t start iteration step %d:"%iter
212 eps_e = eps_e_old + deps-(dt/2)*i_eta*deviatoric(sigma)
213
214 I1=trace(eps_e)
215 sqI2=length(eps_e)
216 xi=safeDiv(I1,sqI2)
217 i_xi=safeDiv(sqI2,I1)
218 # update damage parameter:
219 m=wherePositive(xi-XI_0)
220 a=sqI2**2*(xi-XI_0)*(m*C_D + (1-m)* C_1)
221 b=(1-m)*(1./C_2)
222
223 alpha=solveODE(alpha_old, a,b, dt)
224 alpha_dot=(alpha-alpha_old)/dt
225 i_eta = clip(2*C_V*alpha_dot,minval=0.)
226
227 if inf(alpha) < -EPSILON*10:
228 raise ValueError,"alpha<0"
229 if sup(alpha) > 1:
230 raise ValueError,"alpha > 1"
231 # step size for the next time step:
232
233 gamma=alpha*GAMMA_M
234 lame=LAME_0
235 mu=MU_0*(1-alpha)
236
237 lame_eff=lame-gamma*i_xi
238 mu_eff=mu-gamma*xi/2.
239
240 print "\talpha = [ %e, %e]"%(inf(alpha),sup(alpha))
241 print "\tmu_eff = [ %e, %e]"%(inf(mu_eff),sup(mu_eff))
242 print "\tlame_eff = [ %e, %e]"%(inf(lame_eff),sup(lame_eff))
243 print "\txi = [ %e, %e]"%(inf(xi),sup(xi))
244 print "\tgamma = [ %e, %e]"%(inf(gamma),sup(gamma))
245
246 if inf(mu_eff) < 0:
247 raise ValueError,"mu_eff<0"
248
249 sigma = 2*mu_eff*eps_e+lame_eff*trace(eps_e)*k3
250
251 if (UPDATE_OPERATOR) :
252 pde.setValue(A = mu_eff * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame_eff*k3Xk3)
253 else:
254 pde.setValue(A = mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3)
255
256 pde.setValue(X=-sigma, y=SIGMA_N*dom.getNormal(), r=dt*v0-du)
257
258
259 ddu=pde.getSolution()
260 deps+=symmetric(grad(ddu))
261 du+=ddu
262 norm_ddu=Lsup(ddu)
263 norm_du=Lsup(du)
264 print "\t displacement change update = %e of %e"%(norm_ddu, norm_du)
265 iter+=1
266 print "deps =", inf(deps),sup(deps)
267 n+=1
268 t+=dt
269 #=========== this is a test for triaxial test ===========================
270 print "\tYYY t = ", t
271 a =(SIGMA_N-lame_eff*VMAX*t)/(lame_eff+mu_eff)/2
272 #=========== this is a test for triaxial test ===========================
273 print "\tYYY a = ", meanValue(a)
274 print "\tYYY eps00 = ",meanValue( eps_e[0,0])
275 print "\tYYY eps11 = ",meanValue( eps_e[1,1])
276 print "\tYYY eps22 = num/exact", meanValue(eps_e[2,2]), VMAX*t
277 print "\tYYY eps_kk = num/exact", meanValue(trace(eps_e)), meanValue(VMAX*t+2*a)
278 print "\tYYY sigma11 = num/exact", meanValue(sigma[1,1]), meanValue(lame_eff*(VMAX*t+2*a)+2*mu_eff*a)
279 print "\tYYY sigma22 = num/exact", meanValue(sigma[2,2]), meanValue(lame_eff*(VMAX*t+2*a)+2*mu_eff*VMAX*t)
280 print "\tYYY linear Elastic equivalent num/exact=",meanValue(sigma[2,2]-sigma[0,0]-(sigma_old[2,2]-sigma_old[0,0]))/meanValue(eps_e[2,2]-eps_e_old[2,2]), meanValue(mu_eff*(3*lame_eff+2*mu_eff)/(lame_eff+mu_eff))
281 diagnose.write(("%e,"*9+"\n")%(t, meanValue(-eps_e[2,2]),
282 meanValue(eps_e[1,1]),
283 meanValue(sigma[0,0]-sigma[2,2]),
284 meanValue(mu_eff),
285 meanValue(lame_eff),
286 meanValue(xi),
287 meanValue(gamma),
288 meanValue(alpha),
289 meanValue(alpha_dot)))
290 print "time step %s (t=%s) completed."%(n,t)
291 #
292 # .... visualization
293 #
294 if t>=t_vis or n>n_vis:
295 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),alpha=alpha, I1=trace(eps_e), I2=length(eps_e)**2, xi=safeDiv(trace(eps_e),length(eps_e)))
296 print "visualization file %d for time step %e generated."%(counter_vis,t)
297 counter_vis+=1
298 t_vis+=DT_VIS
299 n_vis+=DN_VIS
300 #
301 # control time step size:
302 #
303 ss=sup(length(deps))
304 if ss>0:
305 dt_new=DEPS_MAX/ss*dt
306 print "\ttime step size to control strain increment %s."%(dt_new,)
307 else:
308 dt_new=dt
309 if dt_old != None:
310 dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
311 norm_alpha=Lsup(alpha)
312 fac=Lsup(dd_alpha)
313 if norm_alpha > 0: fac*=1./norm_alpha
314 if fac>0:
315 error=fac*0.5*dt**2
316 print "\testimated local error for time step size %e is %e"%(dt,error)
317 dt_new=min(dt_new,sqrt(ODE_TOL*2/fac) )
318 else:
319 dt_new=DT_MAX
320
321 dt_new=min(max(dt_new,dt/5),dt*5,DT_MAX) # aviod rapit changes
322 print "\tINFO: new time step size %e"%dt_new
323 dt, dt_old=dt_new, dt
324 # dom.setX(dom.getX()+du)

  ViewVC Help
Powered by ViewVC 1.1.26