/[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 821 by gross, Mon Aug 28 08:36:30 2006 UTC revision 825 by gross, Tue Aug 29 06:59:27 2006 UTC
# Line 96  class Mechanics(Model): Line 96  class Mechanics(Model):
96            k3=kronecker(self.domain)            k3=kronecker(self.domain)
97            # set new thermal stress increment            # set new thermal stress increment
98            if self.temperature:            if self.temperature:
99               self.dthermal_stress=self.bulk_modulus*self.self.expansion_coefficient*(self.temperature-self.__temperature_old)               self.deps_therm=self.self.expansion_coefficient*(self.temperature-self.__temperature_old)
100            else:            else:
101               self.dthermal_stress=0.               self.deps_therm=0.
102            # set PDE coefficients:            # set PDE coefficients:
103            self.__pde.setValue(A=self.S)            self.__pde.setValue(A=self.S)
104            self.__pde.setValue(X=self.stress_old-self.dthermal_stress*k3)            self.__pde.checkSymmetry()
105              self.__pde.setValue(X=self.stress_old-self.bulk_modulus*self.deps_therm*k3)
106            if self.internal_force: self.__pde.setValue(Y=self.internal_force)            if self.internal_force: self.__pde.setValue(Y=self.internal_force)
107            if self.external_force: self.__pde.setValue(y=self.external_force)            if self.external_force: self.__pde.setValue(y=self.external_force)
108            self.__pde.setValue(r=self.prescribed_velocity, \            self.__pde.setValue(r=self.prescribed_velocity, \
# Line 131  class Mechanics(Model): Line 132  class Mechanics(Model):
132             """             """
133             accept all the values:             accept all the values:
134             """             """
135             self.displacement=self.__displacement_old             self.__displacement_old=self.displacement
136             self.stress=self.stress_old             self.stress_old=self.stress
137             self.velocity=self.__velocity_old             self.__velocity_old=self.velocity
138             self.temperature=self.__temperature_old             self.__temperature_old=self.temperature
139    
140        def getSafeTimeStepSize(self,dt):        def getSafeTimeStepSize(self,dt):
141             """             """
# Line 161  class DruckerPrager(Mechanics): Line 162  class DruckerPrager(Mechanics):
162             self.declareParameter(plastic_stress=0.,             self.declareParameter(plastic_stress=0.,
163                                   friction_parameter=0.,                                   friction_parameter=0.,
164                                   dilatancy_parameter=0.,                                   dilatancy_parameter=0.,
165                                   shear_length=1.e15)                                   shear_length=1.e15,
166                                     hardening=0.)
167    
168        def doInitialization(self):        def doInitialization(self):
169             """             """
# Line 170  class DruckerPrager(Mechanics): Line 172  class DruckerPrager(Mechanics):
172             super(DruckerPrager, self).doInitialization()             super(DruckerPrager, self).doInitialization()
173             self.__plastic_stress_old=self.plastic_stress             self.__plastic_stress_old=self.plastic_stress
174             self.__tau_y_old=self.shear_length             self.__tau_y_old=self.shear_length
175               self.__hardening_old=self.hardening
176    
177        def doStepPreprocessing(self,dt):        def doStepPreprocessing(self,dt):
178              """              """
# Line 180  class DruckerPrager(Mechanics): Line 183  class DruckerPrager(Mechanics):
183              """              """
184              super(DruckerPrager, self).doStepPreprocessing(dt)              super(DruckerPrager, self).doStepPreprocessing(dt)
185              self.plastic_stress=self.__plastic_stress_old              self.plastic_stress=self.__plastic_stress_old
186                self.shear_length=self.__tau_y_old
187                self.hardening=self.__hardening_old
188    
189        def doStep(self,dt):        def doStep(self,dt):
190             G=self.shear_modulus             G=self.shear_modulus
# Line 187  class DruckerPrager(Mechanics): Line 192  class DruckerPrager(Mechanics):
192             alpha=self.friction_parameter             alpha=self.friction_parameter
193             beta=self.dilatancy_parameter             beta=self.dilatancy_parameter
194             tau_Y=self.shear_length             tau_Y=self.shear_length
195             if self.__plastic_stress_old:             h=self.hardening
               dps=self.plastic_stress-self.__plastic_stress_old  
               h=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))  
            else:  
               h=0  
196             # set new tangential operator:             # set new tangential operator:
197             self.S=self.getTangentialTensor(self.stress,             self.S=self.getTangentialTensor(self.stress,
198                                             tau_Y,G,K,alpha,beta,h)                                             tau_Y,G,K,alpha,beta,h)
# Line 201  class DruckerPrager(Mechanics): Line 202  class DruckerPrager(Mechanics):
202             # update stresses:             # update stresses:
203             self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,             self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,
204                                                          self.velocity*dt,                                                          self.velocity*dt,
205                                                          self.dthermal_stress,tau_Y,G,K,alpha,beta,h)                                                          self.deps_therm,tau_Y,G,K,alpha,beta,h)
206               if self.__plastic_stress_old:
207                  dps=self.plastic_stress-self.__plastic_stress_old
208                  self.hardening=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
209               else:
210                  self.hardening=0
211    
212        def doStepPostprocessing(self,dt):        def doStepPostprocessing(self,dt):
213            super(DruckerPrager, self).doStepPostprocessing(dt)            super(DruckerPrager, self).doStepPostprocessing(dt)
214            self.plastic_stress=self.__plastic_stress_old=self.plastic_stress            self.__plastic_stress_old=self.plastic_stress
215              self.__tau_y_old=self.shear_length
216              self.__hardening_old=self.hardening
217    
218        def getNewStress(self,s,gamma_p,du,ds_therm,tau_Y,G,K,alpha,beta,h):        def getNewStress(self,s,gamma_p,du,deps_therm,tau_Y,G,K,alpha,beta,h):
219              k3=kronecker(self.domain)              k3=kronecker(self.domain)
220              dt=1.              dt=1.
221              g=grad(du)              g=grad(du)
222              D=symmetric(g)              D=symmetric(g)
223              W=nonsymmetric(g)              W=nonsymmetric(g)
224              s_e=s+ds_therm+dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \              s_e=s+K*deps_therm*k3 +dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \
225                                 +2*nonsymmetric(matrix_mult(W,s)))                                       +2*nonsymmetric(matrix_mult(W,s)))
226              p_e=-1./3*trace(s_e)              p_e=-1./3*trace(s_e)
227              s_e_dev=s_e+p_e*k3              s_e_dev=s_e+p_e*k3
228              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 235  class DruckerPrager(Mechanics): Line 243  class DruckerPrager(Mechanics):
243             chi=whereNonNegative(tau-alpha*p-tau_Y)             chi=whereNonNegative(tau-alpha*p-tau_Y)
244             sXk3=outer(s,k3)             sXk3=outer(s,k3)
245             k3Xk3=outer(k3,k3)             k3Xk3=outer(k3,k3)
246             tmp=G*s_dev/tau             tmp=G*s_dev/(tau+tau_Y*1.e-15*whereZero(tau,1.e-15))
247             S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \             S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
248                 + (K-2./3*G)*k3Xk3 \                 + (K-2./3*G)*k3Xk3 \
249                 + sXk3-swap_axes(sXk3,1,3) \                 + sXk3-swap_axes(sXk3,1,3) \
250                 + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \                 + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \
251                        -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2))                        -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2))
252                 # - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\                 - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\
253               # print S
254             return S             return S
Legend:
Removed from v.821  
changed lines
  Added in v.825
  ViewVC Help
Powered by ViewVC 1.1.26