/[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 828 by gross, Tue Aug 29 09:12:33 2006 UTC revision 829 by gross, Thu Aug 31 01:31:15 2006 UTC
# Line 18  class Mechanics(Model): Line 18  class Mechanics(Model):
18        @ivar internal_force: =Data()        @ivar internal_force: =Data()
19        @ivar external_force: =Data()        @ivar external_force: =Data()
20        @ivar prescribed_velocity: =Data()        @ivar prescribed_velocity: =Data()
21        @ivar location_prescribed_velocity: =Data()        @ivar location_fixed_velocity: =Data()
22        @ivar temperature:  = None        @ivar temperature:  = None
23        @ivar expansion_coefficient:  = 0.        @ivar expansion_coefficient:  = 0.
24        @ivar bulk_modulus: =1.        @ivar bulk_modulus: =1.
# Line 44  class Mechanics(Model): Line 44  class Mechanics(Model):
44                                 velocity=None, \                                 velocity=None, \
45                                 internal_force=Data(), \                                 internal_force=Data(), \
46                                 external_force=Data(), \                                 external_force=Data(), \
47                                 prescribed_velocity=Data(), \                                 location_fixed_velocity=Data(), \
                                location_prescribed_velocity=Data(), \  
48                                 temperature = None, \                                 temperature = None, \
49                                 expansion_coefficient = 0., \                                 expansion_coefficient = 0., \
50                                 bulk_modulus=2., \                                 bulk_modulus=2., \
# Line 60  class Mechanics(Model): Line 59  class Mechanics(Model):
59             if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))             if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
60             if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))             if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
61             if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))             if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
62             self.__pde=LinearPDE(self.domain)             # open PDE:
            self.__pde.setSymmetryOn()  
            self.__displacement_old=self.displacement  
63             self.stress_old=self.stress             self.stress_old=self.stress
64             self.__velocity_old=self.velocity             self.__velocity_old=self.velocity
65             self.__temperature_old=self.temperature             self.__temperature_old=self.temperature
66                         self.__displacement_old=self.displacement
67               # get node cooridnates:
68               self.__x=self.domain.getX()
69               self.domain.setX(self.__x+self.displacement)
70               # open PDE:
71               self.__pde=LinearPDE(self.domain)
72               self.__pde.setSymmetryOn()
73        def doStepPreprocessing(self,dt):        def doStepPreprocessing(self,dt):
74              """              """
75              step up pressure iteration              step up pressure iteration
# Line 76  class Mechanics(Model): Line 79  class Mechanics(Model):
79              """              """
80              self.__iter=0              self.__iter=0
81              self.__diff=self.UNDEF_DT              self.__diff=self.UNDEF_DT
82              # set new values:              # set initial values for the iteration:
             self.displacement=self.__displacement_old  
             self.stress=self.stress_old  
83              self.velocity=self.__velocity_old              self.velocity=self.__velocity_old
84              self.temperature=self.__temperature_old              self.temperature=self.__temperature_old
85              self.__velocity_last=self.velocity              self.displacement=self.__displacement_old+self.velocity*dt
86    
87                # set temperature strain:
88                if self.temperature:
89                   self.deps_therm=self.self.expansion_coefficient*(self.temperature-self.__temperature_old)
90                else:
91                   self.deps_therm=0.
92    
93                # update geometry
94                self.domain.setX(self.__x+self.displacement)
95                self.stress=self.stress_old
96    
97        def doStep(self,dt):        def doStep(self,dt):
98            """            """
# Line 102  class Mechanics(Model): Line 113  class Mechanics(Model):
113            # set PDE coefficients:            # set PDE coefficients:
114            self.__pde.setValue(A=self.S)            self.__pde.setValue(A=self.S)
115            self.__pde.checkSymmetry()            self.__pde.checkSymmetry()
116            self.__pde.setValue(X=self.stress_old-self.bulk_modulus*self.deps_therm*k3)            self.__pde.setValue(X=-self.stress+self.stress_old-self.bulk_modulus*self.deps_therm*k3)
117            if self.internal_force: self.__pde.setValue(Y=self.internal_force)            if self.internal_force: self.__pde.setValue(Y=self.internal_force)
118            if self.external_force: self.__pde.setValue(y=self.external_force)            if self.external_force: self.__pde.setValue(y=self.external_force)
119            self.__pde.setValue(r=self.prescribed_velocity, \            self.__pde.setValue(q=self.location_fixed_velocity)
                               q=self.location_prescribed_velocity)  
120            # solve the PDE:            # solve the PDE:
121            self.__pde.setTolerance(self.rel_tol/100.)            self.__pde.setTolerance(self.rel_tol**2)
122            self.velocity=self.__pde.getSolution(verbose=True)            du=self.__pde.getSolution(verbose=True)
123            # calculate convergence indicators:            self.velocity+=du/dt
124            self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff            # update geometry
           self.__velocity_last=self.velocity  
125            self.displacement=self.__displacement_old+dt*self.velocity            self.displacement=self.__displacement_old+dt*self.velocity
126              self.domain.setX(self.__x+self.displacement)
127    
128            self.__iter+=1            self.__iter+=1
129            if self.debug:            if self.debug:
130               for i in range(self.domain.getDim()):               for i in range(self.domain.getDim()):
131                  self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))                  self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
132               self.trace("velocity increment %s"%self.__diff)               for i in range(self.domain.getDim()):
133            if self.__iter>2 and diff_old<self.__diff:                  self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
               raise IterationDivergenceError,"no improvement in stress iteration"  
           if self.__iter>self.max_iter:  
               raise IterationDivergenceError,"Maximum number of iterations steps reached"  
134    
135        def terminateIteration(self):        def terminateIteration(self):
136            """iteration is terminateIterationd if relative pressure change is less then rel_tol"""            """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
137            return self.__diff<=self.rel_tol*Lsup(self.velocity)+self.abs_tol            if self.__iter>2 and self.diff_last<self.diff:
138                  raise IterationDivergenceError,"no improvement in stress iteration"
139              if self.__iter>self.max_iter:
140                  raise IterationDivergenceError,"Maximum number of iterations steps reached"
141              return self.diff<=self.rel_tol*1.e-2*Lsup(self.stress)+self.abs_tol
142    
143        def doStepPostprocessing(self,dt):        def doStepPostprocessing(self,dt):
144             """             """
145             accept all the values:             accept all the values:
146             """             """
147             self.__displacement_old=self.displacement             self.__displacement_old=self.displacement
148             self.stress_old=self.stress             self.__velocity_old,self.__velocity_old=self.velocity,self.__velocity_old
            self.__velocity_old=self.velocity  
149             self.__temperature_old=self.temperature             self.__temperature_old=self.temperature
150               self.stress_old=self.stress
151    
152        def getSafeTimeStepSize(self,dt):        def getSafeTimeStepSize(self,dt):
153             """             """
154             returns new step size             returns new step size
155             """             """
156             d=Lsup(self.velocity-self.__velocity_old)/dt             d=Lsup(self.velocity-self.__velocity_old)*Lsup(self.displacement)
157             if d>0:             if d>0:
158                 return Lsup(self.displacement)/d                return dt/d*self.rel_tol
159             else:             else:
160                 return self.UNDEF_DT                return self.UNDEF_DT
161    
162    
163    
# Line 162  class DruckerPrager(Mechanics): Line 174  class DruckerPrager(Mechanics):
174             self.declareParameter(plastic_stress=0.,             self.declareParameter(plastic_stress=0.,
175                                   friction_parameter=0.,                                   friction_parameter=0.,
176                                   dilatancy_parameter=0.,                                   dilatancy_parameter=0.,
177                                   shear_length=1.e15,                                   shear_length=1.e15)
                                  hardening=0.)  
   
178        def doInitialization(self):        def doInitialization(self):
179             """            """
180             initialize model            """
181             """            super(DruckerPrager, self).doInitialization()
            super(DruckerPrager, self).doInitialization()  
            self.__plastic_stress_old=self.plastic_stress  
            self.__tau_y_old=self.shear_length  
            self.__hardening_old=self.hardening  
182    
183        def doStepPreprocessing(self,dt):        def doStepPreprocessing(self):
184              """            """
185              step up pressure iteration            """
186              super(DruckerPrager, self).doStepPreprocessing(dt)
187              # set initial guess for iteration:
188              self.shear_length=self.__shear_length_old
189              self.plastic_stress=self.__plastic_stress_old
190              # set initial stress:
191              self.updateStress()
192    
193              if run within a time dependend problem extrapolation of pressure from previous time steps is used to  
             get an initial guess (that needs some work!!!!!!!)  
             """  
             super(DruckerPrager, self).doStepPreprocessing(dt)  
             self.plastic_stress=self.__plastic_stress_old  
             self.shear_length=self.__tau_y_old  
             self.hardening=self.__hardening_old  
194    
195        def doStep(self,dt):  #================================
196              
197              # set new tangential operator:
198              self.S=self.getTangentialTensor()
199              # do the update step:
200    
201              super(DruckerPrager, self).doStep(dt)
202              # update stresses:
203              self.updateStress()
204    
205          def doStepPostprocessing(self,dt):
206              super(DruckerPrager, self).doStepPostprocessing(dt)
207              self.__plastic_stress_old=self.plastic_stress
208              self.__shear_length_old=self.shear_length
209    
210          def setStress(self):
211             G=self.shear_modulus             G=self.shear_modulus
212             K=self.bulk_modulus             K=self.bulk_modulus
213             alpha=self.friction_parameter             alpha=self.friction_parameter
214             beta=self.dilatancy_parameter             beta=self.dilatancy_parameter
            tau_Y=self.shear_length  
215             h=self.hardening             h=self.hardening
216             # set new tangential operator:             k3=kronecker(self.domain)
217             self.S=self.getTangentialTensor(self.stress,             # elastic trial stress:
218                                             tau_Y,G,K,alpha,beta,h)             g=grad(self.velocity)
219             # do the update step:             D=symmetric(g)
220             super(DruckerPrager, self).doStep(dt)             W=nonsymmetric(g)
221               s_e=self.stress_old+K*self.deps_therm*k3+ \
222             # update stresses:                        dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \
223             self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,                            +2*symmetric(matrix_mult(W,self.stress_old)))
224                                                          self.velocity*dt,             p_e=-1./3*trace(s_e)
225                                                          self.deps_therm,tau_Y,G,K,alpha,beta,h)             s_e_dev=s_e+p_e*k3
226             if self.__plastic_stress_old:             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
227                dps=self.plastic_stress-self.__plastic_stress_old             # yield conditon for elastic trial stress:
228                self.hardening=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))             F=tau_e-alpha*p_e-self.shear_length
229             else:             self.__chi=whereNonNegative(F)
230                self.hardening=0             # plastic stress increment:
231               l=self.__chi*F/(h+G+beta*K)
232        def doStepPostprocessing(self,dt):             self.__tau=tau_e-G*l
233            super(DruckerPrager, self).doStepPostprocessing(dt)             self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev+(p_e+l*beta*K)*k3
234            self.__plastic_stress_old=self.plastic_stress             self.plastic_stress=self.plastic_stress_old+l
235            self.__tau_y_old=self.shear_length             # update hardening
236            self.__hardening_old=self.hardening             self.hardening=(self.shear_length-self.__shear_length_old)/(l+self.abs_tol*whereZero(l))
237    
238        def getNewStress(self,s,gamma_p,du,deps_therm,tau_Y,G,K,alpha,beta,h):        def setTangentialTensor(self):
239              k3=kronecker(self.domain)             G=self.shear_modulus
240              dt=1.             K=self.bulk_modulus
241              g=grad(du)             alpha=self.friction_parameter
242              D=symmetric(g)             beta=self.dilatancy_parameter
243              W=nonsymmetric(g)             tau_Y=self.shear_length
244              s_e=s+K*deps_therm*k3 +dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \             chi=self.__chi
245                                       +2*nonsymmetric(matrix_mult(W,s)))             h=self.hardening
246              p_e=-1./3*trace(s_e)             tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))
             s_e_dev=s_e+p_e*k3  
             tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))  
             F=tau_e-alpha*p_e-tau_Y  
             chi=whereNonNegative(F)  
             l=chi*F/(h+G+beta*K)  
             s=(1.-l*G/tau_e)*s_e_dev+(p_e+l*beta*K)*k3  
             gamma_p=gamma_p+l  
             return s, gamma_p  
               
247    
       def getTangentialTensor(self,s,tau_Y,G,K,alpha,beta,h):  
            d=self.domain.getDim()  
248             k3=kronecker(Function(self.domain))             k3=kronecker(Function(self.domain))
249             p=-1./d*trace(s)             s_dev=self.stress+trace(self.stress)*(k3/d)
250             s_dev=s+p*k3             sXk3=outer(self.stress_old,k3)
            tau=sqrt(1./2*inner(s_dev,s_dev))  
            chi=whereNonNegative(tau-alpha*p-tau_Y)  
            sXk3=outer(s,k3)  
251             k3Xk3=outer(k3,k3)             k3Xk3=outer(k3,k3)
            tmp=G*s_dev/(tau+tau_Y*1.e-15*whereZero(tau,1.e-15))  
252             S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \             S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
253                 + (K-2./3*G)*k3Xk3 \                 + (K-2./3*G)*k3Xk3 \
254                 + sXk3-swap_axes(sXk3,1,3) \                 + (sXk3-swap_axes(sXk3,1,3)) \
255                 + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \                 + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \
256                        -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2)) \                        -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2)) \
257                 - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\                 - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)
            # print S  
258             return S             return S
Legend:
Removed from v.828  
changed lines
  Added in v.829
  ViewVC Help
Powered by ViewVC 1.1.26