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

  ViewVC Help
Powered by ViewVC 1.1.26