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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26