modellib/py_src/mechanics.py

# $Id:$

__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""

from esys.escript import *
from esys.escript.modelframe import Model,IterationDivergenceError
from esys.escript.linearPDEs import LinearPDE

class Mechanics(Model):
      """
      base class for mechanics models in updated lagrangean framework

      @ivar domain: domain (in)
      @ivar internal_force: =Data()
      @ivar external_force: =Data()
      @ivar prescribed_velocity: =Data()
      @ivar location_prescribed_velocity: =Data()
      @ivar temperature:  = None
      @ivar expansion_coefficient:  = 0.
      @ivar bulk_modulus: =1.
      @ivar shear_modulus: =1.
      @ivar rel_tol: =1.e-3
      @ivar abs_tol: =1.e-15
      @ivar max_iter: =10
      @ivar displacement: =None
      @ivar stress: =None
      @ivar velocity: =None
      """
      def __init__(self,debug=False):
         """
         set up the model
         
         @param debug: debug flag
         @type debug: C{bool}
         """
         super(Mechanics, self).__init__(self,debug=debug)
         self.declareParameter(domain=None, \
                               displacement=None, \
                               stress=None, \
                               velocity=None, \
                               internal_force=Data(), \
                               external_force=Data(), \
                               prescribed_velocity=Data(), \
                               location_prescribed_velocity=Data(), \
                               temperature = None, \
                               expansion_coefficient = 0., \
                               bulk_modulus=2., \
                               shear_modulus=1., \
                               rel_tol=1.e-3,abs_tol=1.e-15,max_iter=10)
         self.__iter=0

      def doInitialization(self):
           """
           initialize model
           """
           if not self.displacement: self.displacement=Vector(0.,ContinuousFunction(self.domain))
           if not self.velocity: self.velocity=Vector(0.,ContinuousFunction(self.domain))
           if not self.stress: self.stress=Tensor(0.,ContinuousFunction(self.domain))
           self.__pde=LinearPDE(self.domain)
           self.__pde.setSymmetryOn()
           self.__displacement_old=self.displacement
           self.stress_old=self.stress
           self.__velocity_old=self.velocity
           self.__temperature_old=self.temperature
           
      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            self.__iter=0
            self.__diff=self.UNDEF_DT
            # set new values:
            self.displacement=self.__displacement_old
            self.stress=self.stress_old
            self.velocity=self.__velocity_old
            self.temperature=self.__temperature_old
            self.__velocity_last=self.velocity

      def doStep(self,dt):
          """

          performs an iteration step of the penalty method.
          IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
     
          requires self.S to be set 
          updates the thermal stress increment

          """
          k3=kronecker(self.domain)
          # set new thermal stress increment
          if self.temperature:
             self.deps_therm=self.self.expansion_coefficient*(self.temperature-self.__temperature_old)
          else:
             self.deps_therm=0.
          # set PDE coefficients:
          self.__pde.setValue(A=self.S)
          self.__pde.checkSymmetry()
          self.__pde.setValue(X=self.stress_old-self.bulk_modulus*self.deps_therm*k3)
          if self.internal_force: self.__pde.setValue(Y=self.internal_force)
          if self.external_force: self.__pde.setValue(y=self.external_force)
          self.__pde.setValue(r=self.prescribed_velocity, \
                              q=self.location_prescribed_velocity)
          # solve the PDE:
          self.__pde.setTolerance(self.rel_tol/100.)
          self.velocity=self.__pde.getSolution(verbose=True)
          # calculate convergence indicators:
          self.__diff,diff_old=Lsup(self.velocity-self.__velocity_last),self.__diff
          self.__velocity_last=self.velocity
          self.displacement=self.__displacement_old+dt*self.velocity
          self.__iter+=1
          if self.debug:
             for i in range(self.domain.getDim()):
                self.trace("velocity %d range %e:%e"%(i,inf(self.velocity[i]),sup(self.velocity[i])))
             self.trace("velocity increment %s"%self.__diff)
          if self.__iter>2 and diff_old<self.__diff:
              raise IterationDivergenceError,"no improvement in stress iteration"
          if self.__iter>self.max_iter:
              raise IterationDivergenceError,"Maximum number of iterations steps reached"

      def terminateIteration(self):
          """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
          return self.__diff<=self.rel_tol*Lsup(self.velocity)+self.abs_tol

      def doStepPostprocessing(self,dt):
           """
           accept all the values:
           """
           self.__displacement_old=self.displacement
           self.stress_old=self.stress
           self.__velocity_old=self.velocity
           self.__temperature_old=self.temperature

      def getSafeTimeStepSize(self,dt):
           """
           returns new step size
           """
           d=Lsup(self.velocity-self.__velocity_old)/dt
           if d>0:
               return Lsup(self.displacement)/d
           else:
               return self.UNDEF_DT


class DruckerPrager(Mechanics):
      """

      """

      def __init__(self,debug=False):
           """
           set up model
           """
           super(DruckerPrager, self).__init__(debug=debug)
           self.declareParameter(plastic_stress=0.,
                                 friction_parameter=0.,
                                 dilatancy_parameter=0.,
                                 shear_length=1.e15,
                                 hardening=0.)

      def doInitialization(self):
           """
           initialize model
           """
           super(DruckerPrager, self).doInitialization()
           self.__plastic_stress_old=self.plastic_stress
           self.__tau_y_old=self.shear_length
           self.__hardening_old=self.hardening

      def doStepPreprocessing(self,dt):
            """
            step up pressure iteration

            if run within a time dependend problem extrapolation of pressure from previous time steps is used to
            get an initial guess (that needs some work!!!!!!!)
            """
            super(DruckerPrager, self).doStepPreprocessing(dt)
            self.plastic_stress=self.__plastic_stress_old
            self.shear_length=self.__tau_y_old
            self.hardening=self.__hardening_old

      def doStep(self,dt):
           G=self.shear_modulus
           K=self.bulk_modulus
           alpha=self.friction_parameter
           beta=self.dilatancy_parameter
           tau_Y=self.shear_length
           h=self.hardening
           # set new tangential operator:
           self.S=self.getTangentialTensor(self.stress,
                                           tau_Y,G,K,alpha,beta,h)
           # do the update step:
           super(DruckerPrager, self).doStep(dt)

           # update stresses:
           self.stress,self.plastic_stress=self.getNewStress(self.stress_old,self.__plastic_stress_old,
                                                        self.velocity*dt,
                                                        self.deps_therm,tau_Y,G,K,alpha,beta,h)
           if self.__plastic_stress_old:
              dps=self.plastic_stress-self.__plastic_stress_old
              self.hardening=(tau_Y-self.__tau_y_old)/(dps+self.abs_tol*whereZero(dps))
           else:
              self.hardening=0

      def doStepPostprocessing(self,dt):
          super(DruckerPrager, self).doStepPostprocessing(dt)
          self.__plastic_stress_old=self.plastic_stress
          self.__tau_y_old=self.shear_length
          self.__hardening_old=self.hardening

      def getNewStress(self,s,gamma_p,du,deps_therm,tau_Y,G,K,alpha,beta,h):
            k3=kronecker(self.domain)
            dt=1.
            g=grad(du)
            D=symmetric(g)
            W=nonsymmetric(g)
            s_e=s+K*deps_therm*k3 +dt*(2*G*D+(K-2./3*G)*trace(D)*k3 \
                                     +2*nonsymmetric(matrix_mult(W,s)))
            p_e=-1./3*trace(s_e)
            s_e_dev=s_e+p_e*k3
            tau_e=sqrt(1./2*inner(s_e_dev,s_e_dev))
            F=tau_e-alpha*p_e-tau_Y
            chi=whereNonNegative(F)
            l=chi*F/(h+G+beta*K)
            s=(1.-l*G/tau_e)*s_e_dev+(p_e+l*beta*K)*k3
            gamma_p=gamma_p+l
            return s, gamma_p
            

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