/[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 3275 - (show annotations)
Thu Oct 14 06:15:27 2010 UTC (8 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 10310 byte(s)
As it turns out the tangential operator implementation is not a good idea.
This a check in to save the last version of this approach before starting using a different method.


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/450.
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
105
106 #diagnose=FileWriter("diagnose.csv",append=False)
107 diagnose=FileWriter("/mnt/D/Cache/diagnose.csv",append=False)
108 #===================================
109 S=0.5*XI_0*((2.*MU_0+3.*LAME_0)/(3.-XI_0**2) + LAME_0)
110 GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
111 def solveODE(u0, a, b, dt):
112 """
113 solves du/dt=a*exp(b*u) u(t=0)=u0 and return approximation at t=dt
114
115 we sue backwards Euler by solving u-u0 = dt * a*exp(b*u)
116 with newton scheme u <- u - (u-u0 - dt * a*exp(b*u)) / (1-dt*a*b*exp(b*u))
117 """
118 u=u0.copy()
119 norm_du=1.
120 norm_u=0.
121 n=0
122 while norm_du > ODE_ITER_TOL * norm_u:
123 H=-dt*a*exp(b*u)
124 du=-(u-u0+H)/(1+b*H)
125 u+=du
126 norm_du = Lsup(du)
127 norm_u = Lsup(u)
128 n+=1
129 if n>ODE_ITER_MAX: raise ValueError,"ODE iteration failed."
130 print "\tODE iteration completed after %d steps."%(n,)
131 return u
132
133 #======================
134 t=0 # time stamp
135 n=0 # time step counter
136 dt=DT # current time step size
137 t_vis=0
138 n_vis=0
139 counter_vis=0
140 mkDir(VIS_DIR)
141 #=========================
142 #
143 # set up domain
144 #
145 if DIM==2:
146 dom=Rectangle(NE_L,NE_H,l0=L,l1=H,order=-1,optimize=True)
147 else:
148 dom=Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H,order=-1,optimize=True)
149
150 BBOX=boundingBox(dom)
151 DIM=dom.getDim()
152 x=dom.getX()
153 #
154 # initial values:
155 #
156 sigma=Tensor(0.,Function(dom))
157 eps_e=Tensor(0.,Function(dom))
158 if CASE==2 or CASE==3:
159 alpha=ALPHA_0*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
160 else:
161 alpha=Scalar(ALPHA_0,Function(dom))
162 u=Vector(0.,ContinuousFunction(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 u0=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 u0[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 u0[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 #
187 # let the show begin:
188 #
189 k3=kronecker(DIM)
190 k3Xk3=outer(k3,k3)
191 alpha_old=alpha
192 dt_old=None
193 if CASE == 1:
194 diagnose.write("t, -e22, s00-s22, alpha\n")
195 else:
196 diagnose.write("t, eps00, eps11, error0, sigma11/tau\n")
197 while t<T_END:
198
199 print "start time step %d"%(n+1,)
200
201 I1=trace(eps_e)
202 sqI2=length(eps_e)
203 xi=safeDiv(I1,sqI2)
204 i_xi=safeDiv(sqI2,I1)
205 # update damage parameter:
206 m=wherePositive(xi-XI_0)
207 a=sqI2**2*(xi-XI_0)*(m*C_D + (1-m)* C_1)
208 b=(1-m)*(1./C_2)
209
210 alpha, alpha_old, alpha_oold =solveODE(alpha, a,b, dt), alpha, alpha_old
211 alpha_dot=(alpha-alpha_old)/dt
212
213 if inf(alpha) < -EPSILON*10:
214 raise ValueError,"alpha<0"
215 if sup(alpha) > 1:
216 raise ValueError,"alpha > 1"
217 # step size for the next time step:
218
219 gamma=alpha*GAMMA_M
220 lame=LAME_0
221 mu=MU_0*(1-alpha)
222
223 lame_eff=lame-gamma*i_xi
224 mu_eff=mu-gamma*xi
225 print "\talpha = [ %e, %e]"%(inf(alpha),sup(alpha))
226 print "\tmu_eff = [ %e, %e]"%(inf(mu_eff),sup(mu_eff))
227 print "\tlame_eff = [ %e, %e]"%(inf(lame_eff),sup(lame_eff))
228 print "\txi = [ %e, %e]"%(inf(xi),sup(xi))
229 print "\tgamma = [ %e, %e]"%(inf(gamma),sup(gamma))
230
231 if inf(mu_eff) < 0:
232 raise ValueError,"mu_eff<0"
233
234 if CASE == 1 or CASE ==2:
235 if t>0:
236 diagnose.write(("%e,"*4+"\n")%(t, meanValue(-eps_e[2,2]), meanValue(sigma[0,0]-sigma[2,2]), meanValue(alpha)))
237 else:
238 if n>1:
239 print "\t eps00 = ", eps_e[0,0]
240 print "\t eps10 = ", eps_e[1,0]
241 print "\t eps01 = ", eps_e[0,1]
242 print "\t eps11 = ", eps_e[1,1]
243 print "\t eps11/eps00 = ", eps_e[1,1]/eps_e[0,0]
244 print" error:",Lsup(((2*mu-gamma*safeDiv(I1,sqI2))*eps_e[1,1]-(gamma*sqI2-lame*I1)) )/Lsup(length(sigma))
245 print" error:",Lsup(sigma[1,1])/Lsup(length(sigma))
246 diagnose.write("%e, %e, %e, %e, %e\n"%(t,inf(eps_e[0,0]),inf(eps_e[1,1]), Lsup(((2*mu-gamma*safeDiv(I1,sqI2))*eps_e[1,1]-(gamma*sqI2-lame*I1)) )/Lsup(length(sigma)), Lsup(sigma[1,1])/Lsup(length(sigma))))
247
248
249 i_eta = clip(2*C_V*alpha_dot,minval=0.)
250
251 # update strain:
252
253 H = safeDiv(I1,sqI2**2)*eps_e-k3
254 S= -GAMMA_M*I1*i_xi * k3 -(MU_0+GAMMA_M*xi/2) * eps_e
255 pde.setValue(A = mu_eff * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame_eff*k3Xk3 + gamma*i_xi * outer(H,H),
256 X=- (sigma + dt * (S * alpha_dot + mu_eff * i_eta * gamma * safeDiv(length(deviatoric(eps_e))**2,sqI2) * H ) ) ,
257 r=dt*u0, y=SIGMA_N*dom.getNormal())
258
259 du=pde.getSolution()
260 deps=symmetric(grad(du))
261 print "\tYYY deps =", inf(deps),sup(deps)
262
263 eps_e, eps_e_old=eps_e+deps-dt/2*i_eta*deviatoric(sigma), eps_e
264 sigma, sigma_old=2*mu_eff*eps_e+lame_eff*trace(eps_e)*k3, sigma
265 print "time step %s (t=%s) completed."%(n,t)
266 n+=1
267 t+=dt
268 print "\tYYY t = ", t
269 a =(SIGMA_N-lame_eff*VMAX*t)/(lame_eff+mu_eff)/2
270 print "\tYYY a = ", inf(a)
271 print "\tYYY eps00 = ",inf( eps_e[0,0])
272 print "\tYYY eps11 = ",inf( eps_e[1,1])
273 print "\tYYY eps22 = ", inf(eps_e[2,2]), VMAX*t
274 print "\tYYY trace = ", inf(trace(eps_e)), inf(VMAX*t+2*a)
275 print "\tYYY sigma11 = ", inf(sigma[1,1]), inf(lame_eff*(VMAX*t+2*a)+2*mu_eff*a)
276 print "\tYYY sigma22 = ", inf(sigma[2,2]), inf(lame_eff*(VMAX*t+2*a)+2*mu_eff*VMAX*t)
277 print "\tYYY linear Elastic equivalent =",inf(sigma[2,2]-sigma[0,0]-(sigma_old[2,2]-sigma_old[0,0]))/inf(eps_e[2,2]-eps_e_old[2,2]), inf(mu_eff*(3*lame_eff+2*mu_eff)/(lame_eff+mu_eff))
278 #
279 # .... visualization
280 #
281 if t>=t_vis or n>n_vis:
282 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)))
283 print "visualization file %d for time step %e generated."%(counter_vis,t)
284 counter_vis+=1
285 t_vis+=DT_VIS
286 n_vis+=DN_VIS
287 #
288 # control time step size:
289 #
290 ss=sup(length(deps))
291 if ss>0:
292 dt_new=DEPS_MAX/ss*dt
293 print "\ttime step size to control strain increment %s."%(dt_new,)
294 else:
295 dt_new=dt
296 if dt_old != None:
297 dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
298 norm_alpha=Lsup(alpha)
299 fac=Lsup(dd_alpha)
300 if norm_alpha > 0: fac*=1./norm_alpha
301 if fac>0:
302 error=fac*0.5*dt**2
303 print "\testimated local error for time step size %e is %e"%(dt,error)
304 dt_new=min(dt_new,sqrt(ODE_TOL*2/fac) )
305 else:
306 dt_new=DT_MAX
307
308 dt_new=min(max(dt_new,dt/5),dt*5,DT_MAX) # avid rapit changes
309 print "\tINFO: new time step size %e"%dt_new
310 dt, dt_old=dt_new, dt
311 # dom.setX(dom.getX()+du)

  ViewVC Help
Powered by ViewVC 1.1.26