# Diff of /trunk/escript/py_src/rheologies.py

revision 2299 by caltinay, Wed Dec 17 03:08:58 2008 UTC revision 2300 by gross, Wed Mar 11 08:17:57 2009 UTC
# Line 37  import util Line 37  import util
37  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE
38  from pdetools import Defect, NewtonGMRES, ArithmeticTuple  from pdetools import Defect, NewtonGMRES, ArithmeticTuple
39
40    class PowerLaw(object):
41        """
42        this implements the power law for a composition of a set of materials where the viscosity eta of each material is given by a
43        power law relationship of the form
44
45        M{eta=eta_N*(tau/tau_t)**(1.-1./power)}
46
47        where tau is equivalent stress and eta_N, tau_t and power are given constant. Moreover an elastic component can be considered.
48        Moreover tau meets the Drucker-Prager type yield condition
49
50        M{tau <= tau_Y + friction * pressure}
51
52        where gamma_dot is the equivalent.
53        """
54        def __init__(self, numMaterials=1,verbose=False):
55             """
56             initializes a power law
57
58             @param numMaterials: number of materials
59             @type numMaterials: C{int}
60             @param verbose: if C{True} some informations are printed.
61             @type verbose: C{bool}
62             """
63             if numMaterials<1:
64                raise ValueError,"at least one material must be defined."
65             self.__numMaterials=numMaterials
66             self.__eta_N=[None for i in xrange(self.__numMaterials)]
67             self.__tau_t=[1. for i in xrange(self.__numMaterials)]
68             self.__power=[1. for i in xrange(self.__numMaterials)]
69             self.__tau_Y=None
70             self.__friction=0
71             self.__mu=None
72             self.__verbose=verbose
73             self.setEtaTolerance()
74        #===========================================================================
75        def getNumMaterials(self):
76             """
77             returns the numebr of materials
78
79             @return: number of materials
80             @rtype: C{int}
81             """
82             return self.__numMaterials
83        def validMaterialId(self,id=0):
84             """
85             checks if a given material id is valid
86
87             @param id: a material id
88             @type id: C{int}
89             @return: C{True} is the id is valid
90             @rtype: C{bool}
91             """
92             return 0<=id and id<self.getNumMaterials()
93        def setEtaTolerance(self,rtol=util.sqrt(util.EPSILON)):
94             """
95             sets the relative tolerance for the effectice viscosity.
96
97             @param rtol: relative tolerance
98             @type rtol: positive C{float}
99             """
100             if rtol<=0:
101                 raise ValueError,"rtol needs to positive."
102             self.__rtol=rtol
103        def getEtaTolerance(self):
104             """
105             returns the relative tolerance for the effectice viscosity.
106
107             @return: relative tolerance
108             @rtype rtol: positive C{float}
109             """
110             return self.__rtol
111        #===========================================================================
112        def setDruckerPragerLaw(self,tau_Y=None,friction=0):
113              """
114              Sets the parameters for the Drucker-Prager model.
115
116              @param tau_Y: yield stress
117              @param friction: friction coefficient
118              """
119              self.__tau_Y=tau_Y
120              self.__friction=friction
121        def getFriction(self):
122             """
123             returns the friction coefficient
124
125             @return: friction coefficient
126             """
127             return self.__friction
128        def getTauY(self):
129             """
130             returns the yield stress
131
132             @return: the yield stress
133             """
134             return self.__tau_Y
135        #===========================================================================
136        def getElasticShearModulus(self):
137            """
138            returns the elastic shear modulus.
139
140            @return: elastic shear modulus
141            """
142            return self.__mu
143        def setElasticShearModulus(self,mu=None):
144            """
145            Sets the elastic shear modulus.
146
147            @param mu: elastic shear modulus
148            """
149            self.__mu=mu
150        #===========================================================================
151        def getPower(self, id=None):
152             """
153             returns the power in the power law
154
155             @param id:  if present, the power for material C{id} is returned.
156             @type id: C{int}
157             @return: the list of the powers for all matrials is returned. If C{id} is present only the power for material C{id} is returned.
158             """
159             if id == None:
160                return self.__power
161             else:
162                if self.validMaterialId(id):
163                  return self.__power[id]
164                else:
165                  raise ValueError,"Illegal material id %s."%id
166        def getEtaN(self, id=None):
167             """
168             returns the viscosity
169
170             @param id:  if present, the viscosity for material C{id} is returned.
171             @type id: C{int}
172             @return: the list of the viscosities for all matrials is returned. If C{id} is present only the viscosity for material C{id} is returned.
173             """
174             if id == None:
175                return self.__eta_N
176             else:
177                if self.validMaterialId(id):
178                  return self.__eta_N[id]
179                else:
180                 raise ValueError,"Illegal material id %s."%id
181        def getTauT(self, id=None):
182             """
183             returns the transition stress
184
185             @param id:  if present, the transition stress for material C{id} is returned.
186             @type id: C{int}
187             @return: the list of the transition stresses for all matrials is returned. If C{id} is present only the transition stress for material C{id} is returned.
188             """
189             if id == None:
190                return self.__tau_t
191             else:
192                if self.validMaterialId(id):
193                  return self.__tau_t[id]
194                else:
195                  raise ValueError,"Illegal material id %s."%id
196
197        def setPowerLaw(self,eta_N, id=0, tau_t=1, power=1):
198              """
199              Sets the power-law parameters for material id
200
201              @param id: material id
202              @type id: C{int}
203              @param eta_N: viscosity for tau=tau_t
204              @param tau_t: transition stress
205              @param power: power law coefficient
206              """
207              if self.validMaterialId(id):
208                 self.__eta_N[id]=eta_N
209                 self.__power[id]=power
210                 self.__tau_t[id]=tau_t
211              else:
212                  raise ValueError,"Illegal material id %s."%id
213
214        def setPowerLaws(self,eta_N, tau_t, power):
215              """
216              Sets the parameters of the power-law for all materials.
217
218              @param eta_N: list of viscosities for tau=tau_t
219              @param tau_t: list of transition stresses
220              @param power: list of power law coefficient
221              """
222              if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
223                  raise ValueError,"%s materials are expected."%self.__numMaterials
224              for i in xrange(self.__numMaterials):
225                   self.setPowerLaw(id=i, eta_N=eta_N[i],tau_t=tau_t[i],power=power[i])
226
227        #===========================================================================
228        def getEtaEff(self,gamma_dot, eta0=None, pressure=0,dt=None, iter_max=10):
229             """
230             returns the effective viscosity eta_eff such that
231
232             M{tau=eta_eff * gamma_dot}
233
234             by solving a non-linear problem for tau.
235
236             @param gamma_dot: equivalent strain gamma_dot
237             @param eta0: initial guess for the effective viscosity (e.g from a previous time step). If not present, an initial guess is calculated.
238             @param pressure: pressure used to calculate yield condition
239             @param dt: time step size. only needed if elastic component is considered.
240             @type dt: positive C{float} if present
241             @param iter_max: maximum number of iteration steps.
242             @type iter_max: C{int}
243             @return: effective viscosity.
244             """
245             numMaterial=self.getNumMaterials()
246             s=[1.-1./p for p in self.getPower() ]
247             eta_N=self.getEtaN()
248             tau_t=self.getTauT()
249             mu=self.getElasticShearModulus()
250             fric=self.getFriction()
251             tau_Y=self.getTauY()
252             if eta0==None:
253                 theta=0.
254                 for i in xrange(numMaterial):
255                      inv_eta_i=0**s[i]/eta_N[i]
256                      theta=theta+inv_eta_i
257                 if util.inf(theta)<=0:
258                     raise ValueError,"unable to set positive initial guess for eta_eff. Most likely no power law with power 1 set."
259                 eta_eff=1./theta
260             else:
261                 if util.inf(eta0)<=0:
262                     raise ValueError,"initial guess for eta_eff is not positive."
263                 eta_eff=eta0
264
265             if mu !=None and dt == None:
266                 raise ValueError,"Time stepsize dt must be given."
267             if dt !=None:
268                 if dt<=0: raise ValueEror,"time step size must be positive."
269             if tau_Y==None:
270                 eta_max=None
271             else:
272                  if util.inf(pressure)*fric<0:
273                       raise ValueError,"pressure needs to be non-negative."
274                  eta_max=(tau_Y+fric*pressure)/(gamma_dot+util.whereZero(gamma_dot)/util.DBLE_MAX)
275             rtol=self.getEtaTolerance()
276             iter =0
277             converged=False
278             tau=eta_eff*gamma_dot
279             if self.__verbose: print "Start calculation of eta_eff (tolerance = %s)\ninitial max eta_eff = %s, tau = %s."%(rtol,util.Lsup(eta_eff),util.Lsup(tau))
280             while not converged:
281                 if iter>max(iter_max,1):
282                    raise RunTimeError,"tolerance not reached after %s steps."%max(iter_max,1)
283                 #===========================================
284                 theta=0. # =1/eta
285                 omega=0. # = tau*theta'= eta'*tau/eta**2
286                 if mu !=None: theta=1./(dt*mu)
287                 for i in xrange(numMaterial):
288                      inv_eta_i=(tau/tau_t[i])**s[i]/eta_N[i]
289                      theta=theta+inv_eta_i
290                      omega=omega+s[i]*inv_eta_i
291                 #===========================================
292                 eta_eff, eta_eff_old=util.clip(eta_eff*(theta+omega)/(eta_eff*theta**2+omega),maxval=eta_max), eta_eff
293                 tau=eta_eff*gamma_dot
294                 d=util.Lsup(eta_eff-eta_eff_old)
295                 l=util.Lsup(eta_eff)
296                 iter+=1
297                 if self.__verbose: print "step %s: correction = %s, max eta_eff = %s, max tau= %s"%(iter, d, l,util.Lsup(tau))
298                 converged= d<= rtol* l
299             return eta_eff
300
301    #====================================================================================================================================
302
303  class IncompressibleIsotropicKelvinFlow(Defect):  class IncompressibleIsotropicKelvinFlow(Defect):
304        """        """
305        This class implements the rheology of an isotropic Kelvin material.        This class implements the rheology of an isotropic Kelvin material.
# Line 68  class IncompressibleIsotropicKelvinFlow( Line 331  class IncompressibleIsotropicKelvinFlow(
331           @param useJaumannStress: C{True} if Jaumann stress is used           @param useJaumannStress: C{True} if Jaumann stress is used
332                                    (not supported yet)                                    (not supported yet)
333           """           """
if numMaterials<1:
raise ValueError,"at least one material must be defined."
334           super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)           super(IncompressibleIsotropicKelvinFlow, self).__init__(**kwargs)
335           self.__domain=domain           self.__domain=domain
336           self.__t=t           self.__t=t
337           self.__vol=util.integrate(1.,Function(self.__domain))           self.__vol=util.integrate(1.,Function(self.__domain))
338           self.__useJaumannStress=useJaumannStress           self.__useJaumannStress=useJaumannStress
self.__numMaterials=numMaterials
self.__eta_N=[None for i in xrange(self.__numMaterials)]
self.__tau_t=[None for i in xrange(self.__numMaterials)]
self.__power=[1 for i in xrange(self.__numMaterials)]
self.__tau_Y=None
self.__friction=None
self.__mu=None
339           self.__v_boundary=Vector(0,Solution(self.__domain))           self.__v_boundary=Vector(0,Solution(self.__domain))
340           #=======================           #=======================
341           # state variables:           # state variables:
# Line 160  class IncompressibleIsotropicKelvinFlow( Line 414  class IncompressibleIsotropicKelvinFlow(
414            """            """
415            return self.__t            return self.__t
416
def setPowerLaw(self,id,eta_N, tau_t=None, power=1):
"""
Sets the power-law parameters for material q.
"""
if id<0 or id>=self.__numMaterials:
raise ValueError,"Illegal material id %s."%id
self.__eta_N[id]=eta_N
self.__power[id]=power
self.__tau_t[id]=tau_t

def setPowerLaws(self,eta_N, tau_t, power):
"""
Sets the parameters of the power-law for all materials.
"""
if len(eta_N)!=self.__numMaterials or len(tau_t)!=self.__numMaterials or len(power)!=self.__numMaterials:
raise ValueError,"%s materials are expected."%self.__numMaterials
for i in xrange(self.__numMaterials):
self.setPowerLaw(i,eta_N[i],tau_t[i],power[i])

def setDruckerPragerLaw(self,tau_Y=None,friction=0):
"""
Sets the parameters for the Drucker-Prager model.
"""
self.__tau_Y=tau_Y
self.__friction=friction

def setElasticShearModulus(self,mu=None):
"""
Sets the elastic shear modulus.
"""
self.__mu=mu
417        def setExternals(self, F=None, f=None, q=None, v_boundary=None):        def setExternals(self, F=None, f=None, q=None, v_boundary=None):
418            """            """
419            Sets externals.            Sets externals.

Legend:
 Removed from v.2299 changed lines Added in v.2300