# Contents of /trunk/modellib/py_src/mechanics.py

Revision 3981 - (show annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years ago) by jfenwick
File MIME type: text/x-python
File size: 12003 byte(s)
```First pass of updating copyright notices
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2012 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 16 __copyright__="""Copyright (c) 2003-2012 by University of Queensland 17 http://www.uq.edu.au 18 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 __url__= 22 23 #from esys.escript import * 24 from esys.escript.modelframe import Model,IterationDivergenceError 25 from esys.escript.linearPDEs import LinearPDE 26 27 class Mechanics(Model): 28 """ 29 base class for mechanics models in updated lagrangean framework 30 31 :ivar domain: domain (in) 32 :ivar internal_force: =Data() 33 :ivar external_force: =Data() 34 :ivar prescribed_velocity: =Data() 35 :ivar location_prescribed_velocity: =Data() 36 :ivar temperature: = None 37 :ivar expansion_coefficient: = 0. 38 :ivar bulk_modulus: =1. 39 :ivar shear_modulus: =1. 40 :ivar rel_tol: =1.e-3 41 :ivar abs_tol: =1.e-15 42 :ivar max_iter: =10 43 :ivar displacement: =None 44 :ivar stress: =None 45 """ 46 SAFTY_FACTOR_ITERATION=1./100. 47 def __init__(self,**kwargs): 48 """ 49 set up the model 50 51 :keyword debug: debug flag 52 :type debug: ``bool`` 53 """ 54 super(Mechanics, self).__init__(self,**kwargs) 55 self.declareParameter(domain=None, \ 56 displacement=None, \ 57 stress=None, \ 58 velocity=None, \ 59 internal_force=None, \ 60 external_force=None, \ 61 prescribed_velocity=None, \ 62 location_prescribed_velocity=None, \ 63 temperature = None, \ 64 expansion_coefficient = 0., \ 65 bulk_modulus=2., \ 66 shear_modulus=1., \ 67 rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10) 68 self.__iter=0 69 70 def doInitialization(self): 71 """ 72 initialize model 73 """ 74 if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain)) 75 if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain)) 76 if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain)) 77 if not self.internal_force: self.internal_force = Data() 78 if not self.external_force: self.external_force = Data() 79 if not self.prescribed_velocity: self.prescribed_velocity = Data() 80 if not self.location_prescribed_velocity: self.location_prescribed_velocity =Data() 81 # save the old values: 82 self.__stress_safe=self.stress 83 self.__temperature_safe=self.temperature 84 self.__displacement_safe=self.displacement 85 self.__velocity_safe=self.velocity 86 self.__velocity_old=None 87 self.__old_dt=None 88 self.__very_old_dt=None 89 # get node cooridnates and apply initial displacement 90 self.__x=self.domain.getX() 91 self.domain.setX(self.__x+self.displacement) 92 # open PDE: 93 self.__pde=LinearPDE(self.domain) 94 self.__pde.setSolverMethod(self.__pde.DIRECT) 95 self.__solver_options=self.__pde.getSolverOptions() 96 self.__solver_options.setSolverMethod(self.__solver_options.DIRECT) 97 self.__solver_options.setVerbosity(self.debug) 98 99 # self.__pde.setSymmetryOn() 100 101 def doStepPreprocessing(self,dt): 102 """ 103 step up pressure iteration 104 105 if run within a time dependend problem extrapolation of pressure from previous time steps is used to 106 get an initial guess (that needs some work!!!!!!!) 107 """ 108 # reset iteration counters: 109 self.__iter=0 110 self.__diff=self.UNDEF_DT 111 # set initial guesses for the iteration: 112 self.displacement=self.__displacement_safe 113 self.stress=self.__stress_safe 114 self.velocity=self.__velocity_safe 115 # update geometry 116 self.domain.setX(self.__x+self.displacement) 117 118 def doStep(self,dt): 119 """ 120 """ 121 self.__iter+=1 122 k3=kronecker(self.domain) 123 # set new thermal stress increment 124 if self.temperature == None: 125 self.deps_th=0. 126 else: 127 self.deps_th=self.expansion_coefficient*(self.temperature-self.__temperature_safe) 128 # set PDE coefficients: 129 self.__pde.setValue(A=self.S) 130 self.__pde.setValue(X=-self.stress-self.bulk_modulus*self.deps_th*k3) 131 if self.internal_force: self.__pde.setValue(Y=self.internal_force) 132 if self.external_force: self.__pde.setValue(y=self.external_force) 133 self.__pde.setValue(q=self.location_prescribed_velocity, \ 134 r=Data()) 135 if not self.prescribed_velocity.isEmpty() and self.__iter==1: 136 self.__pde.setValue(r=dt*self.prescribed_velocity) 137 # solve the PDE: 138 self.__solver_options.setTolerance(self.rel_tol**2) 139 self.du=self.__pde.getSolution() 140 # update geometry 141 self.displacement=self.displacement+self.du 142 self.domain.setX(self.__x+self.displacement) 143 self.velocity=(self.displacement-self.__displacement_safe)/dt 144 145 if self.debug: 146 for i in range(self.domain.getDim()): 147 self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i]))) 148 for i in range(self.domain.getDim()): 149 self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i]))) 150 for i in range(self.domain.getDim()): 151 self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i]))) 152 self.__stress_last=self.stress 153 154 def terminateIteration(self): 155 """iteration is terminateIterationd if relative pressure change is less then rel_tol""" 156 if self.__iter>self.max_iter: 157 raise IterationDivergenceError("Maximum number of iterations steps reached") 158 if self.__iter==0: 159 self.__diff=self.UNDEF_DT 160 else: 161 self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff 162 s_sup=Lsup(self.stress) 163 self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff)) 164 if self.__iter>2 and diff_safe0: 183 return 1./a 184 else: 185 return self.UNDEF_DT 186 187 188 189 class DruckerPrager(Mechanics): 190 """ 191 192 """ 193 194 def __init__(self,**kwargs): 195 """ 196 set up model 197 """ 198 super(DruckerPrager, self).__init__(**kwargs) 199 self.declareParameter(plastic_stress=0., 200 hardening=0., 201 friction_parameter=0., 202 dilatancy_parameter=0., 203 shear_length=1.e15) 204 def doInitialization(self): 205 """ 206 """ 207 super(DruckerPrager, self).doInitialization() 208 self.__plastic_stress_safe=self.plastic_stress 209 self.__shear_length_safe=self.shear_length 210 self.__hardening_safe=self.hardening 211 self.__chi_safe=0 212 self.__tau_safe=0 213 214 def doStepPreprocessing(self,dt): 215 """ 216 """ 217 super(DruckerPrager, self).doStepPreprocessing(dt) 218 # set initial guess for iteration: 219 self.shear_length=self.__shear_length_safe 220 self.plastic_stress=self.__plastic_stress_safe 221 self.hardening=self.__hardening_safe 222 self.__chi=self.__chi_safe 223 self.__tau=self.__tau_safe 224 225 def doStep(self,dt): 226 # set new tangential operator: 227 self.setTangentialTensor() 228 # do the update step: 229 super(DruckerPrager, self).doStep(dt) 230 # update stresses: 231 self.setStress() 232 233 def doStepPostprocessing(self,dt): 234 super(DruckerPrager, self).doStepPostprocessing(dt) 235 self.__plastic_stress_safe=self.plastic_stress 236 self.__shear_length_safe=self.shear_length 237 self.__hardening_safe=self.hardening 238 self.__chi_safe=self.__chi 239 self.__tau_safe=self.__tau 240 241 def setStress(self): 242 d=self.domain.getDim() 243 G=self.shear_modulus 244 K=self.bulk_modulus 245 alpha=self.friction_parameter 246 beta=self.dilatancy_parameter 247 h=self.hardening 248 k3=kronecker(self.domain) 249 # elastic trial stress: 250 g=grad(self.du) 251 D=symmetric(g) 252 W=nonsymmetric(g) 253 s_e=self.stress+K*self.deps_th*k3+ \ 254 2*G*D+(K-2./3*G)*trace(D)*k3 \ 255 +2*symmetric(matrix_mult(W,self.stress)) 256 p_e=-1./d*trace(s_e) 257 s_e_dev=s_e+p_e*k3 258 tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev)) 259 # yield conditon for elastic trial stress: 260 F=tau_e-alpha*p_e-self.shear_length 261 self.__chi=whereNonNegative(F+(self.rel_tol*(self.SAFTY_FACTOR_ITERATION)**2)*self.shear_length) 262 # plastic stress increment: 263 l=self.__chi*F/(h+G+beta*K) 264 self.__tau=tau_e-G*l 265 self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev-(p_e+l*beta*K)*k3 266 self.plastic_stress=self.plastic_stress+l 267 # update hardening 268 self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l)) 269 270 def setTangentialTensor(self): 271 d=self.domain.getDim() 272 G=self.shear_modulus 273 K=self.bulk_modulus 274 alpha=self.friction_parameter 275 beta=self.dilatancy_parameter 276 tau_Y=self.shear_length 277 chi=self.__chi 278 tau=self.__tau 279 h=self.hardening 280 k3=kronecker(Function(self.domain)) 281 282 sXk3=outer(self.stress,k3) 283 k3Xk3=outer(k3,k3) 284 s_dev=self.stress-trace(self.stress)*(k3/d) 285 tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol)) 286 287 self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \ 288 + (K-2./3*G)*k3Xk3 \ 289 + (sXk3-swap_axes(swap_axes(sXk3,1,2),2,3)) \ 290 + 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3) \ 291 -swap_axes(swap_axes(sXk3,0,3),2,3) \ 292 -swap_axes(sXk3,1,2) \ 293 +swap_axes(sXk3,1,3) ) \ 294 - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)