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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3274 by gross, Thu Sep 23 06:53:36 2010 UTC revision 3275 by gross, Thu Oct 14 06:15:27 2010 UTC
# Line 31  import time Line 31  import time
31  # ======================= Default Values ==================================================  # ======================= Default Values ==================================================
32  H=1.*U.m                     # height  H=1.*U.m                     # height
33  L=1.*H                           # length  L=1.*H                           # length
34  NE_H=5                           # number of elements in H-direction.  NE_H=1                          # number of elements in H-direction.
35  NE_L=int(ceil(L*NE_H/H))  NE_L=int(ceil(L*NE_H/H))
36    CASE=1
37    
38  #Boundary conditions:  #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).  #   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).
# Line 45  NE_L=int(ceil(L*NE_H/H)) Line 46  NE_L=int(ceil(L*NE_H/H))
46  #   damage and damage rate  #   damage and damage rate
47    
48    
49    T_END=60000000.0*U.sec                       # end time
50    
51  if True:  if CASE==1 or CASE==2:
52     C_V = 1.e-5/(U.Mega*U.Pa)     C_V = 2.e-5/(U.Mega*U.Pa)
53     C_D = 3/U.sec     C_D = 3/U.sec  
54     C_1 = 1e-12/U.sec     C_1 = 1e-12/U.sec  
55     C_2 = 0.03     C_2 = 0.03
56     XI_0 = -0.56     XI_0 = -0.56
57     LAME_0 = 29e9*U.Pa     LAME_0 = 29e9*U.Pa
58     MU_0 = 19e9*U.Pa     MU_0 = 19e9*U.Pa
59     ALPHA_0 = 0.0     if CASE==2:
60          ALPHA_0 = 0.999
61       else:
62          ALPHA_0 = 0.0
63     RHO = 2800*U.kg/U.m**3     RHO = 2800*U.kg/U.m**3
64     G = 10*U.m/U.sec**2     *0     G = 10*U.m/U.sec**2     *0
    DT=1.*U.sec  
65     SIGMA_N=50.*U.Mega*U.Pa     SIGMA_N=50.*U.Mega*U.Pa
66     DIM=3                               DIM=3                          
67     CASE=1     VMAX=-1.*U.m/U.sec/450.
68     VMAX=-0.1*U.m/U.sec     DT_MAX=50.*U.sec
69     DT_MAX=10.*U.sec     DT=DT_MAX/1000000.
70       xc=[L/2,L/2,H/2]
71       WWW=min(H,L)*0.01
72    
73  else:  else:
74     C_V = 3e-11/U.Pa     C_V = 3e-11/U.Pa
75     C_D = 5/U.sec     C_D = 5/U.sec
76     C_1 = 1e-12/U.sec     C_1 = 1e-12/U.sec
77     C_2 = 0.03     C_2 = 0.03
78     XI_0 = -0.8     XI_0 = -0.8
79     LAME_0 = 46e9*U.Pa     LAME_0 = 46e9*U.Pa
80     MU_0 = 30e9*U.Pa     MU_0 = 30e9*U.Pa
81     ALPHA_0 = 0.1     if CASE==3:
82          ALPHA_0 = 0.01
83       else:
84          ALPHA_0 = 0.0
85     RHO = 2800*U.kg/U.m**3     RHO = 2800*U.kg/U.m**3
86     G = 10*U.m/U.sec**2     *0     G = 10*U.m/U.sec**2     *0
87     DT=1.*U.sec/100.     VMAX=-1*U.m/U.sec
88     VMAX=-0.01*U.m/U.sec     DT_MAX=500.*U.sec
89     DT_MAX=10.*U.sec     DT=DT_MAX/100000.
90     SIGMA_N=0     SIGMA_N=0
91     DIM=2                           DIM=2                      
92     CASE=0     xc=[L/2,H/2]
93       WWW=min(H,L)*0.08
94    
 T_END=15.*U.sec                       # end time  
95  VERBOSE=True  VERBOSE=True
96  DT_VIS=T_END/100                # time distane between two visulaization files  DT_VIS=T_END/100                # time distane between two visulaization files
97  DN_VIS=5                        # maximum counter increment between two visulaization files  DN_VIS=1                        # maximum counter increment between two visulaization files
98  VIS_DIR="results"               # name of the director for vis files  VIS_DIR="results"               # name of the director for vis files
99    
100  ODE_TOL=0.01  ODE_TOL=0.01
101  ODE_ITER_TOL=1.e-8  ODE_ITER_TOL=1.e-8
102  ODE_ITER_MAX=15  ODE_ITER_MAX=15
103    DEPS_MAX=0.01
104    
 diagnose=FileWriter("diagnose.csv",append=False)  
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)  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))  GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
# Line 119  def solveODE(u0, a, b, dt): Line 130  def solveODE(u0, a, b, dt):
130      print "\tODE iteration completed after %d steps."%(n,)      print "\tODE iteration completed after %d steps."%(n,)
131      return u      return u
132    
   
133  #======================  #======================
134  t=0         # time stamp  t=0         # time stamp
135  n=0         # time step counter  n=0         # time step counter
# Line 145  x=dom.getX() Line 155  x=dom.getX()
155  #  #
156  sigma=Tensor(0.,Function(dom))  sigma=Tensor(0.,Function(dom))
157  eps_e=Tensor(0.,Function(dom))  eps_e=Tensor(0.,Function(dom))
158  alpha=Scalar(ALPHA_0,Function(dom))  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))  u=Vector(0.,ContinuousFunction(dom))
163    
164  pde=LinearPDESystem(dom)  pde=LinearPDESystem(dom)
# Line 156  pde.getSolverOptions().setSolverMethod(p Line 169  pde.getSolverOptions().setSolverMethod(p
169    
170  fixed_v_mask=Vector(0,Solution(dom))  fixed_v_mask=Vector(0,Solution(dom))
171  u0=Vector(0.,ContinuousFunction(dom))  u0=Vector(0.,ContinuousFunction(dom))
172  if CASE == 1:  if CASE == 1 or CASE==2:
173      for d in range(DIM):      for d in range(DIM):
174         fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)         fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
175         if d == DIM-1:         if d == DIM-1:
# Line 169  else: Line 182  else:
182             fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)             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             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, r=u0, y=SIGMA_N*dom.getNormal())  pde.setValue(Y=-G*RHO*kronecker(DIM)[DIM-1], q=fixed_v_mask)
186  #  #
187  #  let the show begin:  #  let the show begin:
188  #  #
# Line 178  k3Xk3=outer(k3,k3) Line 191  k3Xk3=outer(k3,k3)
191  alpha_old=alpha  alpha_old=alpha
192  dt_old=None  dt_old=None
193  if CASE == 1:  if CASE == 1:
194     diagnose.write("t,s22-s00, s22-p, e00, e11, alpha, alpha_dot\n")     diagnose.write("t, -e22, s00-s22, alpha\n")
195  else:  else:
196     diagnose.write("t, eps00, eps11, error0, sigma11/tau\n")     diagnose.write("t, eps00, eps11, error0, sigma11/tau\n")
197  while t<T_END:  while t<T_END:
# Line 187  while t<T_END: Line 200  while t<T_END:
200    
201      I1=trace(eps_e)      I1=trace(eps_e)
202      sqI2=length(eps_e)      sqI2=length(eps_e)
     print "II2 =",sqI2  
203      xi=safeDiv(I1,sqI2)      xi=safeDiv(I1,sqI2)
204      i_xi=safeDiv(sqI2,I1)      i_xi=safeDiv(sqI2,I1)
205      # update damage parameter:      # update damage parameter:
# Line 198  while t<T_END: Line 210  while t<T_END:
210      alpha, alpha_old, alpha_oold =solveODE(alpha, a,b, dt), alpha, alpha_old      alpha, alpha_old, alpha_oold =solveODE(alpha, a,b, dt), alpha, alpha_old
211      alpha_dot=(alpha-alpha_old)/dt      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:      # step size for the next time step:
218    
     print "\tXXX alpha = [ %e, %e]"%(inf(alpha),sup(alpha))  
     if sup(alpha)  >1:  
         raise ValueError,"damage parameter > 1 detected."  
   
219      gamma=alpha*GAMMA_M      gamma=alpha*GAMMA_M
220      lame=LAME_0      lame=LAME_0
221      mu=MU_0*(1-alpha)      mu=MU_0*(1-alpha)
222    
223      lame_eff=lame-gamma*i_xi      lame_eff=lame-gamma*i_xi
224      mu_eff=mu-gamma*xi      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))      print "\tmu_eff = [ %e, %e]"%(inf(mu_eff),sup(mu_eff))
227      print "\tlame_eff = [ %e, %e]"%(inf(lame_eff),sup(lame_eff))      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:      if CASE == 1 or CASE ==2:
235        diagnose.write(("%e,"*7+"\n")%(t,meanValue(sigma[2,2]-sigma[0,0]), meanValue(deviatoric(sigma)[2,2]), meanValue(eps_e[0,0]), meanValue(eps_e[1,1]),        if t>0:
236                      meanValue(alpha), meanValue(alpha_dot)))           diagnose.write(("%e,"*4+"\n")%(t, meanValue(-eps_e[2,2]), meanValue(sigma[0,0]-sigma[2,2]),  meanValue(alpha)))
237      else:      else:
238         if n>1:         if n>1:
           print "\t xi = [ %e, %e]"%(inf(xi),sup(xi))  
239            print "\t eps00 = ", eps_e[0,0]            print "\t eps00 = ", eps_e[0,0]
240            print "\t eps10 = ", eps_e[1,0]            print "\t eps10 = ", eps_e[1,0]
241            print "\t eps01 = ", eps_e[0,1]            print "\t eps01 = ", eps_e[0,1]
# Line 229  while t<T_END: Line 245  while t<T_END:
245            print" error:",Lsup(sigma[1,1])/Lsup(length(sigma))            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))))            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    
     if inf(mu_eff) <= 0:  
         raise ValueError,"Material failed due to non-positive mu_eff."  
     if inf(lame_eff) <= 0:  
         raise ValueError,"Material failed due to non-positive lame_eff."  
248    
249      i_eta = clip(2*C_V*alpha_dot,minval=0.)      i_eta = clip(2*C_V*alpha_dot,minval=0.)
250    
251      # update strain:      # update strain:
252    
253      H = safeDiv(I1,sqI2**2)*eps_e-k3      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),      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 * mu_eff * i_eta * gamma * safeDiv(length(deviatoric(eps_e))**2,sqI2) * H ) )                   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()*dt      du=pde.getSolution()
260        deps=symmetric(grad(du))
261        print "\tYYY deps =", inf(deps),sup(deps)
262    
263      eps_e=eps_e+symmetric(grad(du))-dt/2*i_eta*deviatoric(sigma)      eps_e, eps_e_old=eps_e+deps-dt/2*i_eta*deviatoric(sigma), eps_e
264      sigma=2*mu_eff*eps_e+lame_eff*trace(eps_e)*k3      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)      print "time step %s (t=%s) completed."%(n,t)
266      n+=1      n+=1
267      t+=dt      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      #  .... visualization
280      #      #
281      if t>=t_vis or n>n_vis:      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)        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)        print "visualization file %d for time step %e generated."%(counter_vis,t)
284        counter_vis+=1        counter_vis+=1
285        t_vis+=DT_VIS        t_vis+=DT_VIS
# Line 261  while t<T_END: Line 287  while t<T_END:
287      #      #
288      #   control time step size:      #   control time step size:
289      #      #
290      if dt_old == None:      ss=sup(length(deps))
291            dt_new=dt      if ss>0:
292           dt_new=DEPS_MAX/ss*dt
293           print "\ttime step size to control strain increment %s."%(dt_new,)
294      else:      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))            dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
298            norm_alpha=Lsup(alpha)            norm_alpha=Lsup(alpha)
299            fac=Lsup(dd_alpha)            fac=Lsup(dd_alpha)
300            if norm_alpha > 0: fac*=1./norm_alpha            if norm_alpha > 0: fac*=1./norm_alpha
301            error=fac*0.5*dt**2            if fac>0:
302            print "\testimated local error for time step size %e is %e"%(dt,error)               error=fac*0.5*dt**2
303            dt_new=sqrt(2./3.*ODE_TOL*2/fac)               print "\testimated local error for time step size %e is %e"%(dt,error)
304            dt_new=min(max(dt_new,dt/5),dt*5,DT_MAX) # avid rapit changes               dt_new=min(dt_new,sqrt(ODE_TOL*2/fac) )
305            print "\tINFO: new time step size %e"%dt_new            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      dt, dt_old=dt_new, dt
311        # dom.setX(dom.getX()+du)

Legend:
Removed from v.3274  
changed lines
  Added in v.3275

  ViewVC Help
Powered by ViewVC 1.1.26