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

  ViewVC Help
Powered by ViewVC 1.1.26