/[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 835 by gross, Fri Sep 1 09:05:16 2006 UTC revision 836 by gross, Mon Sep 4 22:37:25 2006 UTC
# Line 232  class DruckerPrager(Mechanics): Line 232  class DruckerPrager(Mechanics):
232             # elastic trial stress:             # elastic trial stress:
233             g=grad(self.du)             g=grad(self.du)
234             D=symmetric(g)             D=symmetric(g)
            print trace(D)  
235             W=nonsymmetric(g)             W=nonsymmetric(g)
236             s_e=self.stress+K*self.deps_th*k3+ \             s_e=self.stress+K*self.deps_th*k3+ \
237                        2*G*D+(K-2./3*G)*trace(D)*k3 \                        2*G*D+(K-2./3*G)*trace(D)*k3 \
238                        +2*symmetric(matrix_mult(W,self.stress))                        +2*symmetric(matrix_mult(W,self.stress))
239             print s_e             print "sym (e)= ",Lsup(nonsymmetric(s_e))
240             p_e=-1./d*trace(s_e)             p_e=-1./d*trace(s_e)
241             s_e_dev=s_e+p_e*k3             s_e_dev=s_e+p_e*k3
242             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
# Line 249  class DruckerPrager(Mechanics): Line 248  class DruckerPrager(Mechanics):
248             l=self.__chi*F/(h+G+beta*K)             l=self.__chi*F/(h+G+beta*K)
249             self.__tau=tau_e-G*l             self.__tau=tau_e-G*l
250             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
251               print "sym = ",Lsup(nonsymmetric(self.stress))
252             self.plastic_stress=self.plastic_stress+l             self.plastic_stress=self.plastic_stress+l
253             # update hardening             # update hardening
254             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 278  class DruckerPrager(Mechanics): Line 278  class DruckerPrager(Mechanics):
278             k3Xk3=outer(k3,k3)             k3Xk3=outer(k3,k3)
279             s_dev=self.stress-trace(self.stress)*(k3/d)             s_dev=self.stress-trace(self.stress)*(k3/d)
280             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))
281    
282               print "plast grad= ",Lsup(tmp)
283             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)) \
284                   + (K-2./3*G)*k3Xk3 \                   + (K-2./3*G)*k3Xk3 \
285                   + (sXk3-swap_axes(sXk3,1,3)) \                   + sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)   \
286                   + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \                   + 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3)   \
287                        +swap_axes(sXk3,1,3)-swap_axes(sXk3,0,2)) \                          -swap_axes(swap_axes(sXk3,0,3),2,3)   \
288                            -swap_axes(sXk3,1,2)                  \
289                            +swap_axes(sXk3,1,3)                ) \
290                   - 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)
291               print "sym S:",Lsup(nonsymmetric(self.S))

Legend:
Removed from v.835  
changed lines
  Added in v.836

  ViewVC Help
Powered by ViewVC 1.1.26