/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 11034 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26