/[escript]/trunk/modellib/py_src/mechanics.py
ViewVC logotype

Diff of /trunk/modellib/py_src/mechanics.py

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

revision 837 by gross, Mon Sep 4 22:37:25 2006 UTC revision 838 by gross, Tue Sep 5 22:45:55 2006 UTC
# Line 65  class Mechanics(Model): Line 65  class Mechanics(Model):
65             self.__temperature_safe=self.temperature             self.__temperature_safe=self.temperature
66             self.__displacement_safe=self.displacement             self.__displacement_safe=self.displacement
67             self.__velocity_safe=self.velocity             self.__velocity_safe=self.velocity
68             self.__velocity_last=None             self.__velocity_old=None
69               self.__old_dt=None
70               self.__very_old_dt=None
71             # get node cooridnates and apply initial displacement             # get node cooridnates and apply initial displacement
72             self.__x=self.domain.getX()             self.__x=self.domain.getX()
73             self.domain.setX(self.__x+self.displacement)             self.domain.setX(self.__x+self.displacement)
# Line 88  class Mechanics(Model): Line 90  class Mechanics(Model):
90              self.displacement=self.__displacement_safe              self.displacement=self.__displacement_safe
91              self.stress=self.__stress_safe              self.stress=self.__stress_safe
92              self.velocity=self.__velocity_safe              self.velocity=self.__velocity_safe
93              self.__velocity_last=self.__velocity_safe              self.__velocity_old=self.__velocity_safe
94              print "set self.__velocity_last",dt              self.__very_old_dt, self.__old_dt= self.__old_dt, dt
95              # update geometry              # update geometry
96              self.domain.setX(self.__x+self.displacement)              self.domain.setX(self.__x+self.displacement)
97    
# Line 138  class Mechanics(Model): Line 140  class Mechanics(Model):
140            else:            else:
141               self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff               self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
142               s_sup=Lsup(self.stress)               s_sup=Lsup(self.stress)
143                 # self.__diff,diff_safe=Lsup(self.du),self.__diff
144                 # s_sup=Lsup(self.displacement)
145               self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))               self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
146               if self.__iter>2 and diff_safe<self.__diff:               if self.__iter>2 and diff_safe<self.__diff:
147                   raise IterationDivergenceError,"no improvement in stress iteration"                   raise IterationDivergenceError,"no improvement in stress iteration"
# Line 151  class Mechanics(Model): Line 155  class Mechanics(Model):
155             self.__temperature_safe=self.temperature             self.__temperature_safe=self.temperature
156             self.__stress_safe=self.stress             self.__stress_safe=self.stress
157             self.__velocity_safe=self.velocity             self.__velocity_safe=self.velocity
158               self.__old_dt=dt
159    
160        def getSafeTimeStepSize(self,dt):        def getSafeTimeStepSize(self,dt):
161             """             """
162             returns new step size             returns new step size
163             """             """
164             if self.__velocity_last:             print self.__very_old_dt, self.__old_dt, dt
165                a=Lsup(self.velocity-self.__velocity_last)/dt             if self.__very_old_dt:
166                  a=2*Lsup(self.velocity-self.__velocity_old)/(self.__very_old_dt+self.__old_dt)
167                if a>0:                if a>0:
168                   print "new dt:= ",Lsup(self.displacement)/a*self.rel_tol,dt,a                   print "new dt:= ",Lsup(self.velocity)/a*self.rel_tol,dt,self.__old_dt,a
169                   # return Lsup(self.displacement)/a*self.rel_tol                   # return Lsup(self.displacement)/a*self.rel_tol
170                   return self.UNDEF_DT                   return self.UNDEF_DT
171                else:                else:
# Line 236  class DruckerPrager(Mechanics): Line 242  class DruckerPrager(Mechanics):
242             s_e=self.stress+K*self.deps_th*k3+ \             s_e=self.stress+K*self.deps_th*k3+ \
243                        2*G*D+(K-2./3*G)*trace(D)*k3 \                        2*G*D+(K-2./3*G)*trace(D)*k3 \
244                        +2*symmetric(matrix_mult(W,self.stress))                        +2*symmetric(matrix_mult(W,self.stress))
            print "sym (e)= ",Lsup(nonsymmetric(s_e))  
245             p_e=-1./d*trace(s_e)             p_e=-1./d*trace(s_e)
246             s_e_dev=s_e+p_e*k3             s_e_dev=s_e+p_e*k3
247             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
248             # yield conditon for elastic trial stress:             # yield conditon for elastic trial stress:
249             F=tau_e-alpha*p_e-self.shear_length             F=tau_e-alpha*p_e-self.shear_length
250             self.__chi=whereNonNegative(F)             self.__chi=whereNonNegative(F+(self.rel_tol*(self.SAFTY_FACTOR_ITERATION)**2)*self.shear_length)
            print "chi =",Lsup(self.__chi)  
251             # plastic stress increment:             # plastic stress increment:
252             l=self.__chi*F/(h+G+beta*K)             l=self.__chi*F/(h+G+beta*K)
253             self.__tau=tau_e-G*l             self.__tau=tau_e-G*l
254             self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3             self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3
            print "sym = ",Lsup(nonsymmetric(self.stress))  
255             self.plastic_stress=self.plastic_stress+l             self.plastic_stress=self.plastic_stress+l
256             # update hardening             # update hardening
257             self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))             self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))
# Line 265  class DruckerPrager(Mechanics): Line 268  class DruckerPrager(Mechanics):
268             h=self.hardening             h=self.hardening
269             k3=kronecker(Function(self.domain))             k3=kronecker(Function(self.domain))
270    
            ###  
            p_e=-1./d*trace(self.stress)  
            s_e_dev=self.stress+p_e*k3  
            tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))  
            F=tau_e-alpha*p_e-self.shear_length  
            # print F  
            ###  
   
   
271             sXk3=outer(self.stress,k3)             sXk3=outer(self.stress,k3)
272             k3Xk3=outer(k3,k3)             k3Xk3=outer(k3,k3)
273             s_dev=self.stress-trace(self.stress)*(k3/d)             s_dev=self.stress-trace(self.stress)*(k3/d)
274             tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))             tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))
275    
            print "plast grad= ",Lsup(tmp)  
276             self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \             self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
277                   + (K-2./3*G)*k3Xk3 \                   + (K-2./3*G)*k3Xk3 \
278                   + sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)   \                   + sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)   \
# Line 288  class DruckerPrager(Mechanics): Line 281  class DruckerPrager(Mechanics):
281                          -swap_axes(sXk3,1,2)                  \                          -swap_axes(sXk3,1,2)                  \
282                          +swap_axes(sXk3,1,3)                ) \                          +swap_axes(sXk3,1,3)                ) \
283                   - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)                   - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)
            print "sym S:",Lsup(nonsymmetric(self.S))  

Legend:
Removed from v.837  
changed lines
  Added in v.838

  ViewVC Help
Powered by ViewVC 1.1.26