[cig-commits] r4788 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Oct 11 12:02:26 PDT 2006
Author: willic3
Date: 2006-10-11 12:02:26 -0700 (Wed, 11 Oct 2006)
New Revision: 4788
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
Log:
More progress with effective stress formulation.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2006-10-11 18:39:23 UTC (rev 4787)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_6.f 2006-10-11 19:02:26 UTC (rev 4788)
@@ -299,8 +299,8 @@
& rtimdat,rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
& ierr,errstrng)
c
-c... subroutine to compute the current stress for the time-dependent
-c solution.
+c... subroutine to compute the current values for stress, total strain,
+c and viscous strain increment for the time-dependent solution.
c Note that only the upper triangle is used (or available), as
c dmat is assumed to always be symmetric.
c
@@ -311,6 +311,8 @@
include "ndimens.inc"
include "nconsts.inc"
include "rconsts.inc"
+ integer nrpar,nipar
+ parameter(nrpar=10,nipar=0)
c
c... subroutine arguments
c
@@ -327,21 +329,203 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... local constants
+c
+ double precision diag(6)
+ data diag/one,one,one,zero,zero,zero/
+ double precision eng(6)
+ data eng/one,one,one,half,half,half/
+c
+c... local variables
+c
+ integer i
+ integer ipar(nipar)
+ double precision rpar(nrpar)
+ double precision e,pr,anpwr,emhu,rmu,ae,tf
+ double precision emeantdt,smeant,smeantdt
+ double precision sinv1t,sefft,sinv10,seff0
+ double precision epptdt(6),sdevt(6),sdev0(6),sdevtdt(6)
+ double precision sdev,sdevtp,smean0,sdev0
+ double precision eett(6)
+
+c
c... included variable definitions
c
include "rtimdat_def.inc"
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-c... return error code, as this material is not yet defined
+cdebug write(6,*) "Hello from td_strs_6!"
c
- ierr=101
- errstrng="td_strs_6"
+ ierr=0
c
+c... define material properties
+c
+ e=prop(2)
+ pr=prop(3)
+ anpwr=prop(4)
+ emhu=prop(5)
+ rmu=half*e/(one+pr)
+ ae=(one+pr)/e
+ tf=deltp*(one-alfap)
+c
+c... define stress and strain quantities
+c
+ emeantdt=third*(ee(1)+ee(2)+ee(3))
+ smeantdt=e*emean/(one-two*pr)
+ call invar(sdevt,sinv1t,sefft,state)
+ smeant=third*sinv1t
+ call invar(sdev0,sinv10,seff0,state0)
+ smean0=third*sinv10
+ do i=1,nstr
+ epptdt(i)=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
+ end do
+c
+c... define parameters needed by effective stress function
+c
+ rpar(1)=ae
+ rpar(2)=anpwr
+ rpar(3)=emhu
+ rpar(4)=sprod(epptdt,epptdt)+two*ae*sprod(epptdt,sdev0)+
+ & ae*ae*seff0
+ rpar(5)=two*tf*(sprod(epptdt,sdevt)+ae*sprod(sdevt,sdev0))
+ rpar(6)=sefft
+c
+c... call effective stress driver routine
+c
+ call esfcomp_6(sefftdt,rpar,nrpar,ipar,nipar,ierr,errstrng)
+
+
+ call dcopy(nstr,ee,ione,eett,ione)
+ eett(4)=half*eett(4)
+ eett(5)=half*eett(5)
+ eett(6)=half*eett(6)
+
+c
return
end
c
c
+ subroutine esfcomp_6(sefftdt,rpar,nrpar,ipar,nipar,ierr,errstrng)
+c
+c... driver routine to compute the effective stress at the current
+c time step.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "rconsts.inc"
+ double precision brac
+c... This value is arbitrary, but has worked OK in the past.
+ parameter(brac=5000.0d0)
+c
+c... subroutine arguments
+c
+ integer nrpar,nipar,ierr
+ integer ipar(nipar)
+ double precision sefftdt,rpar(nrpar)
+ character errstrng*(*)
+c
+c... included dimension and type statements
+c
+ include "rtimdat_dim.inc"
+ include "rgiter_dim.inc"
+ include "ntimdat_dim.inc"
+c
+c... external routines
+c
+ external esf_6
+c
+c... local variables
+c
+ double precision add,sefflo,seffhi
+c
+c... included variable definitions
+c
+ include "rtimdat_def.inc"
+ include "rgiter_def.inc"
+ include "ntimdat_def.inc"
+c
+cdebug write(6,*) "Hello from esfcomp_6!"
+c
+c******** decide whether to use current stress or elastic stress
+c computed from current strains as initial value. This
+c would just be scalar product of eppdtd/ae.
+ add=brac*stol*
+
+
+
+c
+c
+ subroutine esf_6(seffcur,f,df,rpar,nrpar,ipar,nipar)
+c
+c... subroutine to compute the effective stress function and it's
+c derivative for a given effective stress for a Maxwell power-law
+c viscoelastic material.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "rconsts.inc"
+c
+c... subroutine arguments
+c
+ integer nrpar,nipar
+ integer ipar(nipar)
+ double precision seffcur,f,df,rpar(nrpar)
+c
+c... included dimension and type statements
+c
+ include "rtimdat_dim.inc"
+c
+c... local variables
+c
+ double precision ae,anpwr,emhu,b,c,sefft
+ double precision f1,tf
+ double precision gamtau,sefftau
+c
+c... included variable definitions
+c
+ include "rtimdat_def.inc"
+c
+cdebug write(6,*) "Hello from esf_6!"
+c
+c
+c... handy constants
+c
+ f1=one-alfap
+ tf=deltp*f1
+c
+c... values from parameter list
+c
+ ae=rpar(1)
+ anpwr=rpar(2)
+ emhu=rpar(3)
+ b=rpar(4)
+ c=rpar(5)
+ sefft=rpar(6)
+ d=tf*sefft
+c
+c... additional values dependent on current effective stress
+c
+ sefftau=f1*sefft+alfap*seffcur
+ gamtau=half*(sefftau/emhu)**(anpwr-one)/emhu
+ a=ae+alfap*deltp*gamtau
+c
+c... effective stress function and it's derivative
+c
+ f=a*a*seffcur*seffcur-b+c*gamtau-d*d*gamtau*gamtau
+ df=two*a*a*seffcur+(two*a*alfap*deltp*seffcur*seffcur+
+ & c-two*d*d*gamtau)*
+ & (half*alfap*(anpwr-one)*(sefftau/emhu)**anpwr-two)/(emhu*emhu)
+c
+ return
+ end
+
+c
+c
subroutine td_strs_mat_6(state,dstate,state0,ee,scur,dmat,prop,
& rtimdat,rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
& ierr,errstrng)
More information about the cig-commits
mailing list