/[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 831 - (show annotations)
Thu Aug 31 08:47:40 2006 UTC (12 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 10366 byte(s)
some more work to get the drucker-prager going
1 # $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 @ivar internal_force: =Data()
19 @ivar external_force: =Data()
20 @ivar prescribed_velocity: =Data()
21 @ivar location_prescribed_velocity: =Data()
22 @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 """
32 SAFTY_FACTOR_ITERATION=1./100.
33 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 prescribed_velocity=Data(), \
48 location_prescribed_velocity=Data(), \
49 temperature = None, \
50 expansion_coefficient = 0., \
51 bulk_modulus=2., \
52 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 # save the old values:
64 self.__stress_old=self.stress
65 self.__temperature_old=self.temperature
66 self.__displacement_old=self.displacement
67 # get node cooridnates and apply initial displacement
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.setSolverMethod(self.__pde.DIRECT)
73 self.__pde.setSymmetryOn()
74
75 def doStepPreprocessing(self,dt):
76 """
77 step up pressure iteration
78
79 if run within a time dependend problem extrapolation of pressure from previous time steps is used to
80 get an initial guess (that needs some work!!!!!!!)
81 """
82 # reset iteration counters:
83 self.__iter=0
84 self.__diff=self.UNDEF_DT
85 # set initial values for the iteration:
86 self.temperature=self.__temperature_old
87 self.displacement=self.__displacement_old
88 self.stress=self.__stress_old
89 # update geometry
90 self.domain.setX(self.__x+self.displacement)
91
92 def doStep(self,dt):
93 """
94 """
95 self.__iter+=1
96 k3=kronecker(self.domain)
97 # set new thermal stress increment
98 if self.temperature:
99 self.deps_th=self.self.expansion_coefficient*(self.temperature-self.__temperature_old)
100 else:
101 self.deps_th=0.
102 # set PDE coefficients:
103 self.__pde.setValue(A=self.S)
104 self.__pde.setValue(X=-self.stress-self.bulk_modulus*self.deps_th*k3)
105 if self.internal_force: self.__pde.setValue(Y=self.internal_force)
106 if self.external_force: self.__pde.setValue(y=self.external_force)
107 self.__pde.setValue(q=self.location_prescribed_velocity, \
108 r=Data())
109 if not self.prescribed_velocity.isEmpty() and self.__iter==1:
110 self.__pde.setValue(r=dt*self.prescribed_velocity)
111 # solve the PDE:
112 self.__pde.setTolerance(self.rel_tol**2)
113 self.du=self.__pde.getSolution(verbose=True)
114 # update geometry
115 self.displacement=self.displacement+self.du
116 self.domain.setX(self.__x+self.displacement)
117
118 if self.debug:
119 for i in range(self.domain.getDim()):
120 self.trace("du %d range %e:%e"%(i,inf(self.du[i]),sup(self.du[i])))
121 for i in range(self.domain.getDim()):
122 self.trace("displacement %d range %e:%e"%(i,inf(self.displacement[i]),sup(self.displacement[i])))
123 self.__stress_last=self.stress
124
125 def terminateIteration(self):
126 """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
127 if self.__iter>self.max_iter:
128 raise IterationDivergenceError,"Maximum number of iterations steps reached"
129 print self.__iter
130 if self.__iter==0:
131 self.__diff=self.UNDEF_DT
132 else:
133 self.__diff,diff_old=Lsup(self.stress-self.__stress_last),self.__diff
134 s_sup=Lsup(self.stress)
135 self.trace("stress max and increment :%e, %e"%(s_sup,self.__diff))
136 if diff_old<self.__diff:
137 raise IterationDivergenceError,"no improvement in stress iteration"
138 return self.__diff<=self.rel_tol*self.SAFTY_FACTOR_ITERATION*s_sup+self.abs_tol
139
140 def doStepPostprocessing(self,dt):
141 """
142 accept all the values:
143 """
144 self.__displacement_old=self.displacement
145 self.__temperature_old=self.temperature
146 self.__stress_old=self.stress
147
148 def getSafeTimeStepSize(self,dt):
149 """
150 returns new step size
151 d=Lsup(self.velocity-self.__velocity_old)*Lsup(self.displacement)
152 if d>0:
153 return dt/d*self.rel_tol
154 else:
155 """
156 return self.UNDEF_DT
157
158
159
160 class DruckerPrager(Mechanics):
161 """
162
163 """
164
165 def __init__(self,debug=False):
166 """
167 set up model
168 """
169 super(DruckerPrager, self).__init__(debug=debug)
170 self.declareParameter(plastic_stress=0.,
171 hardening=0.,
172 friction_parameter=0.,
173 dilatancy_parameter=0.,
174 shear_length=1.e15)
175 def doInitialization(self):
176 """
177 """
178 super(DruckerPrager, self).doInitialization()
179 self.__plastic_stress_old=self.plastic_stress
180 self.__shear_length_old=self.shear_length
181 self.__hardening_old=self.hardening
182 self.__chi_old=0
183 self.__tau_old=0
184
185 def doStepPreprocessing(self,dt):
186 """
187 """
188 super(DruckerPrager, self).doStepPreprocessing(dt)
189 # set initial guess for iteration:
190 self.shear_length=self.__shear_length_old
191 self.plastic_stress=self.__plastic_stress_old
192 self.hardening=self.__hardening_old
193 self.__chi=self.__chi_old
194 self.__tau=self.__tau_old
195
196 def doStep(self,dt):
197 # set new tangential operator:
198 self.setTangentialTensor()
199 # do the update step:
200 super(DruckerPrager, self).doStep(dt)
201 # update stresses:
202 self.setStress()
203
204 def doStepPostprocessing(self,dt):
205 super(DruckerPrager, self).doStepPostprocessing(dt)
206 self.__plastic_stress_old=self.plastic_stress
207 self.__shear_length_old=self.shear_length
208 self.__hardening_old=self.hardening
209 self.__chi_old=self.__chi
210 self.__tau_old=self.__tau
211
212 def setStress(self):
213 d=self.domain.getDim()
214 G=self.shear_modulus
215 K=self.bulk_modulus
216 alpha=self.friction_parameter
217 beta=self.dilatancy_parameter
218 h=self.hardening
219 k3=kronecker(self.domain)
220 # elastic trial stress:
221 g=grad(self.du)
222 D=symmetric(g)
223 W=nonsymmetric(g)
224 s_e=self.stress+K*self.deps_th*k3+ \
225 2*G*D+(K-2./3*G)*trace(D)*k3 \
226 +2*symmetric(matrix_mult(W,self.stress))
227 p_e=-1./d*trace(s_e)
228 s_e_dev=s_e+p_e*k3
229 tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
230 # yield conditon for elastic trial stress:
231 F=tau_e-alpha*p_e-self.shear_length
232 self.__chi=whereNonNegative(F)
233 # plastic stress increment:
234 l=self.__chi*F/(h+G+beta*K)
235 self.__tau=tau_e-G*l
236 self.stress=self.__tau/(tau_e+self.abs_tol*whereZero(tau_e,self.abs_tol))*s_e_dev+(p_e+l*beta*K)*k3
237 self.plastic_stress=self.plastic_stress+l
238 # update hardening
239 self.hardening=(self.shear_length-self.__shear_length_old)/(l+self.abs_tol*whereZero(l))
240
241 def setTangentialTensor(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 tau_Y=self.shear_length
248 chi=self.__chi
249 tau=self.__tau
250 h=self.hardening
251
252 k3=kronecker(Function(self.domain))
253 sXk3=outer(self.stress,k3)
254 k3Xk3=outer(k3,k3)
255 s_dev=self.stress-trace(self.stress)*(k3/d)
256 tmp=G*s_dev/(tau+self.abs_tol*whereZero(tau,self.abs_tol))
257 self.S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
258 + (K-2./3*G)*k3Xk3 \
259 + (sXk3-swap_axes(sXk3,1,3)) \
260 + 1./2*(swap_axes(sXk3,0,3)-swap_axes(sXk3,1,2) \
261 -swap_axes(sXk3,1,3)+swap_axes(sXk3,0,2)) \
262 - outer(chi/(h+G+alpha*beta*K)*(tmp+beta*K*k3),tmp+alpha*K*k3)

  ViewVC Help
Powered by ViewVC 1.1.26