/[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 838 - (hide annotations)
Tue Sep 5 22:45:55 2006 UTC (14 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 11578 byte(s)
drucker prager works now but there is still some work on the stepsize control needed.
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 gross 838 self.__velocity_old=None
69     self.__old_dt=None
70     self.__very_old_dt=None
71 gross 831 # get node cooridnates and apply initial displacement
72 gross 829 self.__x=self.domain.getX()
73     self.domain.setX(self.__x+self.displacement)
74     # open PDE:
75     self.__pde=LinearPDE(self.domain)
76 gross 831 self.__pde.setSolverMethod(self.__pde.DIRECT)
77 gross 829 self.__pde.setSymmetryOn()
78 gross 831
79 gross 814 def doStepPreprocessing(self,dt):
80     """
81     step up pressure iteration
82    
83     if run within a time dependend problem extrapolation of pressure from previous time steps is used to
84     get an initial guess (that needs some work!!!!!!!)
85     """
86 gross 831 # reset iteration counters:
87 gross 814 self.__iter=0
88     self.__diff=self.UNDEF_DT
89 gross 834 # set initial guesses for the iteration:
90 gross 832 self.displacement=self.__displacement_safe
91     self.stress=self.__stress_safe
92     self.velocity=self.__velocity_safe
93 gross 838 self.__velocity_old=self.__velocity_safe
94     self.__very_old_dt, self.__old_dt= self.__old_dt, dt
95 gross 829 # update geometry
96     self.domain.setX(self.__x+self.displacement)
97    
98 gross 814 def doStep(self,dt):
99     """
100     """
101 gross 831 self.__iter+=1
102 gross 814 k3=kronecker(self.domain)
103     # set new thermal stress increment
104     if self.temperature:
105 gross 832 self.deps_th=self.self.expansion_coefficient*(self.temperature-self.__temperature_safe)
106 gross 814 else:
107 gross 831 self.deps_th=0.
108 gross 814 # set PDE coefficients:
109     self.__pde.setValue(A=self.S)
110 gross 831 self.__pde.setValue(X=-self.stress-self.bulk_modulus*self.deps_th*k3)
111 gross 814 if self.internal_force: self.__pde.setValue(Y=self.internal_force)
112     if self.external_force: self.__pde.setValue(y=self.external_force)
113 gross 831 self.__pde.setValue(q=self.location_prescribed_velocity, \
114     r=Data())
115     if not self.prescribed_velocity.isEmpty() and self.__iter==1:
116     self.__pde.setValue(r=dt*self.prescribed_velocity)
117 gross 814 # solve the PDE:
118 gross 829 self.__pde.setTolerance(self.rel_tol**2)
119 gross 831 self.du=self.__pde.getSolution(verbose=True)
120 gross 829 # update geometry
121 gross 831 self.displacement=self.displacement+self.du
122 gross 834 self.domain.setX(self.__x+self.displacement)
123 gross 832 self.velocity=(self.displacement-self.__displacement_safe)/dt
124 gross 829
125 gross 821 if self.debug:
126     for i in range(self.domain.getDim()):
127 gross 831 self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
128 gross 829 for i in range(self.domain.getDim()):
129     self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
130 gross 832 for i in range(self.domain.getDim()):
131     self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
132 gross 831 self.__stress_last=self.stress
133 gross 829
134     def terminateIteration(self):
135     """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
136 gross 814 if self.__iter>self.max_iter:
137     raise IterationDivergenceError,"Maximum number of iterations steps reached"
138 gross 831 if self.__iter==0:
139     self.__diff=self.UNDEF_DT
140     else:
141 gross 832 self.__diff,diff_safe=Lsup(self.stress-self.__stress_last),self.__diff
142 gross 831 s_sup=Lsup(self.stress)
143 gross 838 # self.__diff,diff_safe=Lsup(self.du),self.__diff
144     # s_sup=Lsup(self.displacement)
145 gross 831 self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
146 gross 834 if self.__iter>2 and diff_safe<self.__diff:
147 gross 831 raise IterationDivergenceError,"no improvement in stress iteration"
148     return self.__diff<=self.rel_tol*self.SAFTY_FACTOR_ITERATION*s_sup+self.abs_tol
149 gross 814
150     def doStepPostprocessing(self,dt):
151     """
152     accept all the values:
153     """
154 gross 832 self.__displacement_safe=self.displacement
155     self.__temperature_safe=self.temperature
156     self.__stress_safe=self.stress
157     self.__velocity_safe=self.velocity
158 gross 838 self.__old_dt=dt
159 gross 814
160     def getSafeTimeStepSize(self,dt):
161     """
162     returns new step size
163 gross 832 """
164 gross 838 print self.__very_old_dt, self.__old_dt, dt
165     if self.__very_old_dt:
166     a=2*Lsup(self.velocity-self.__velocity_old)/(self.__very_old_dt+self.__old_dt)
167 gross 832 if a>0:
168 gross 838 print "new dt:= ",Lsup(self.velocity)/a*self.rel_tol,dt,self.__old_dt,a
169 gross 832 # return Lsup(self.displacement)/a*self.rel_tol
170     return self.UNDEF_DT
171     else:
172     return self.UNDEF_DT
173 gross 814 else:
174 gross 832 return self.UNDEF_DT
175 gross 814
176    
177    
178     class DruckerPrager(Mechanics):
179     """
180    
181     """
182    
183     def __init__(self,debug=False):
184     """
185     set up model
186     """
187     super(DruckerPrager, self).__init__(debug=debug)
188     self.declareParameter(plastic_stress=0.,
189 gross 831 hardening=0.,
190 gross 814 friction_parameter=0.,
191     dilatancy_parameter=0.,
192 gross 829 shear_length=1.e15)
193 gross 814 def doInitialization(self):
194 gross 829 """
195     """
196     super(DruckerPrager, self).doInitialization()
197 gross 832 self.__plastic_stress_safe=self.plastic_stress
198     self.__shear_length_safe=self.shear_length
199     self.__hardening_safe=self.hardening
200     self.__chi_safe=0
201     self.__tau_safe=0
202 gross 814
203 gross 831 def doStepPreprocessing(self,dt):
204 gross 829 """
205     """
206     super(DruckerPrager, self).doStepPreprocessing(dt)
207     # set initial guess for iteration:
208 gross 832 self.shear_length=self.__shear_length_safe
209     self.plastic_stress=self.__plastic_stress_safe
210     self.hardening=self.__hardening_safe
211     self.__chi=self.__chi_safe
212     self.__tau=self.__tau_safe
213 gross 814
214 gross 831 def doStep(self,dt):
215 gross 829 # set new tangential operator:
216 gross 831 self.setTangentialTensor()
217 gross 829 # do the update step:
218     super(DruckerPrager, self).doStep(dt)
219     # update stresses:
220 gross 831 self.setStress()
221 gross 829
222     def doStepPostprocessing(self,dt):
223     super(DruckerPrager, self).doStepPostprocessing(dt)
224 gross 832 self.__plastic_stress_safe=self.plastic_stress
225     self.__shear_length_safe=self.shear_length
226     self.__hardening_safe=self.hardening
227     self.__chi_safe=self.__chi
228     self.__tau_safe=self.__tau
229 gross 829
230     def setStress(self):
231 gross 831 d=self.domain.getDim()
232 gross 814 G=self.shear_modulus
233     K=self.bulk_modulus
234     alpha=self.friction_parameter
235     beta=self.dilatancy_parameter
236 gross 829 h=self.hardening
237     k3=kronecker(self.domain)
238     # elastic trial stress:
239 gross 831 g=grad(self.du)
240 gross 829 D=symmetric(g)
241     W=nonsymmetric(g)
242 gross 831 s_e=self.stress+K*self.deps_th*k3+ \
243     2*G*D+(K-2./3*G)*trace(D)*k3 \
244     +2*symmetric(matrix_mult(W,self.stress))
245     p_e=-1./d*trace(s_e)
246 gross 829 s_e_dev=s_e+p_e*k3
247     tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
248     # yield conditon for elastic trial stress:
249     F=tau_e-alpha*p_e-self.shear_length
250 gross 838 self.__chi=whereNonNegative(F+(self.rel_tol*(self.SAFTY_FACTOR_ITERATION)**2)*self.shear_length)
251 gross 829 # plastic stress increment:
252     l=self.__chi*F/(h+G+beta*K)
253     self.__tau=tau_e-G*l
254 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
255 gross 831 self.plastic_stress=self.plastic_stress+l
256 gross 829 # update hardening
257 gross 832 self.hardening=(self.shear_length-self.__shear_length_safe)/(l+self.abs_tol*whereZero(l))
258 gross 829
259     def setTangentialTensor(self):
260 gross 831 d=self.domain.getDim()
261 gross 829 G=self.shear_modulus
262     K=self.bulk_modulus
263     alpha=self.friction_parameter
264     beta=self.dilatancy_parameter
265 gross 814 tau_Y=self.shear_length
266 gross 829 chi=self.__chi
267 gross 831 tau=self.__tau
268 gross 825 h=self.hardening
269 gross 832 k3=kronecker(Function(self.domain))
270 gross 814
271 gross 831 sXk3=outer(self.stress,k3)
272 gross 814 k3Xk3=outer(k3,k3)
273 gross 831 s_dev=self.stress-trace(self.stress)*(k3/d)
274     tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))
275 gross 836
276 gross 834 self.S=G*(swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3)) \
277 gross 831 + (K-2./3*G)*k3Xk3 \
278 gross 836 + sXk3-swap_axes(swap_axes(sXk3,1,2),2,3) \
279     + 1./2*(swap_axes(swap_axes(sXk3,0,2),2,3) \
280     -swap_axes(swap_axes(sXk3,0,3),2,3) \
281     -swap_axes(sXk3,1,2) \
282     +swap_axes(sXk3,1,3) ) \
283 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