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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 9 months ago) by ksteube
File MIME type: text/x-python
File size: 11746 byte(s)
Copyright updated in all python files

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

  ViewVC Help
Powered by ViewVC 1.1.26