/[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 819 - (hide 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 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     @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 gross 819 print self.prescribed_velocity
95     print self.location_prescribed_velocity
96 gross 814 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 gross 816 self.__tau_y_old=self.shear_length
159 gross 814
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 gross 816 h=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
179 gross 814 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