/[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 820 - (show annotations)
Mon Aug 28 06:55:36 2006 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 9566 byte(s)
make the modelframe test run
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 @ivar velocity=None
32 """
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=1., \
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 self.__pde=LinearPDE(self.domain)
64 self.__displacement_old=self.displacement
65 self.stress_old=self.stress
66 self.__velocity_old=self.velocity
67 self.__temperature_old=self.temperature
68
69 def doStepPreprocessing(self,dt):
70 """
71 step up pressure iteration
72
73 if run within a time dependend problem extrapolation of pressure from previous time steps is used to
74 get an initial guess (that needs some work!!!!!!!)
75 """
76 self.__iter=0
77 self.__diff=self.UNDEF_DT
78 # set new values:
79 self.displacement=self.__displacement_old
80 self.stress=self.stress_old
81 self.velocity=self.__velocity_old
82 self.temperature=self.__temperature_old
83 self.__velocity_last=self.velocity
84
85 def doStep(self,dt):
86 """
87
88 performs an iteration step of the penalty method.
89 IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
90
91 requires self.S to be set
92 updates the thermal stress increment
93
94 """
95 k3=kronecker(self.domain)
96 # set new thermal stress increment
97 if self.temperature:
98 self.dthermal_stress=self.bulk_modulus*self.self.expansion_coefficient*(self.temperature-self.__temperature_old)
99 else:
100 self.dthermal_stress=0.
101 # set PDE coefficients:
102 self.__pde.setValue(A=self.S)
103 self.__pde.setValue(X=self.stress_old-self.dthermal_stress*k3)
104 if self.internal_force: self.__pde.setValue(Y=self.internal_force)
105 if self.external_force: self.__pde.setValue(y=self.external_force)
106 print self.prescribed_velocity
107 print self.location_prescribed_velocity
108 self.__pde.setValue(r=self.prescribed_velocity, \
109 q=self.location_prescribed_velocity)
110 # solve the PDE:
111 self.__pde.setTolerance(self.rel_tol/100.)
112 self.velocity=self.__pde.getSolution()
113 # calculate convergence indicators:
114 self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff
115 self.__velocity_last=self.velocity
116 self.displacement=self.__displacement_old+dt*self.velocity
117 self.__iter+=1
118 self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
119 if self.__iter>2 and diff_old<self.__diff:
120 raise IterationDivergenceError,"no improvement in stress iteration"
121 if self.__iter>self.max_iter:
122 raise IterationDivergenceError,"Maximum number of iterations steps reached"
123
124 def terminateIteration(self):
125 """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
126 return self.__diff<=self.rel_tol*Lsup(self.velocity)+self.abs_tol
127
128 def doStepPostprocessing(self,dt):
129 """
130 accept all the values:
131 """
132 self.displacement=self.__displacement_old
133 self.stress=self.stress_old
134 self.velocity=self.__velocity_old
135 self.temperature=self.__temperature_old
136
137 def getSafeTimeStepSize(self,dt):
138 """
139 returns new step size
140 """
141 d=Lsup(self.velocity-self.__velocity_old)/dt
142 if d>0:
143 return Lsup(self.displacement)/d
144 else:
145 return self.UNDEF_DT
146
147
148
149 class DruckerPrager(Mechanics):
150 """
151
152 """
153
154 def __init__(self,debug=False):
155 """
156 set up model
157 """
158 super(DruckerPrager, self).__init__(debug=debug)
159 self.declareParameter(plastic_stress=0.,
160 friction_parameter=0.,
161 dilatancy_parameter=0.,
162 shear_length=0.)
163
164 def doInitialization(self):
165 """
166 initialize model
167 """
168 super(DruckerPrager, self).doInitialization()
169 self.__plastic_stress_old=self.plastic_stress
170 self.__tau_y_old=self.shear_length
171
172 def doStepPreprocessing(self,dt):
173 """
174 step up pressure iteration
175
176 if run within a time dependend problem extrapolation of pressure from previous time steps is used to
177 get an initial guess (that needs some work!!!!!!!)
178 """
179 super(DruckerPrager, self).doStepPreprocessing(dt)
180 self.plastic_stress=self.__plastic_stress_old
181
182 def doStep(self,dt):
183 G=self.shear_modulus
184 K=self.bulk_modulus
185 alpha=self.friction_parameter
186 beta=self.dilatancy_parameter
187 tau_Y=self.shear_length
188 if self.__plastic_stress_old:
189 dps=self.plastic_stress-self.__plastic_stress_old
190 h=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
191 else:
192 h=0
193 # set new tangential operator:
194 self.S=self.getTangentialTensor(self.stress,
195 tau_Y,G,K,alpha,beta,h)
196 # do the update step:
197 super(DruckerPrager, self).doStep(dt)
198
199 # update stresses:
200 self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,
201 self.velocity*dt,
202 self.dthermal_stress,tau_Y,G,K,alpha,beta,h)
203
204 def doStepPostprocessing(self,dt):
205 super(DruckerPrager, self).doStepPostprocessing(dt)
206 self.plastic_stress=self.__plastic_stress_old=self.plastic_stress
207
208 def getNewStress(self,s,gamma_p,du,ds_therm,tau_Y,G,K,alpha,beta,h):
209 k3=kronecker(self.domain)
210 dt=1.
211 g=grad(du)
212 D=symmetric(g)
213 W=nonsymmetric(g)
214 s_e=s+ds_therm+dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \
215 +2*nonsymmetric(matrix_mult(W,s)))
216 p_e=-1./3*trace(s_e)
217 s_e_dev=s_e+p_e*k3
218 tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
219 F=tau_e-alpha*p_e-tau_Y
220 chi=whereNonNegative(F)
221 l=chi*F/(h+G+beta*K)
222 s=(1.-l*G/tau_e)*s_e_dev+(p_e+l*beta*K)*k3
223 gamma_p=gamma_p+l
224 return s, gamma_p
225
226
227 def getTangentialTensor(self,s,tau_Y,G,K,alpha,beta,h):
228 d=self.domain.getDim()
229 k3=kronecker(Function(self.domain))
230 p=-1./d*trace(s)
231 s_dev=s+p*k3
232 tau=sqrt(1./2*inner(s_dev,s_dev))
233 chi=whereNonNegative(tau-alpha*p-tau_Y)
234 sXk3=outer(s,k3)
235 k3Xk3=outer(k3,k3)
236 tmp=G*s_dev/tau
237 S=G*(swap_axes(k3Xk3,1,2)+swap_axes(k3Xk3,1,3)) \
238 + (K-2./3*G)*k3Xk3 \
239 + sXk3-swap_axes(sXk3,1,3) \
240 + 1./2*(swap_axes(sXk3,0,3)+swap_axes(sXk3,1,2) \
241 -swap_axes(sXk3,1,3)-swap_axes(sXk3,0,2))
242 # - chi/(h+G+alpha*beta*K)*outer(tmp+beta*K*k3,tmp+alpha*K*k3)\
243 return S

  ViewVC Help
Powered by ViewVC 1.1.26