[cig-commits] r5271 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Nov 15 12:04:40 PST 2006
Author: willic3
Date: 2006-11-15 12:04:40 -0800 (Wed, 15 Nov 2006)
New Revision: 5271
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f
Log:
Initial incomplete set of routines for generalized Maxwell model.
Things should be OK through matinit computation.
I still need to finish the rest of it.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f 2006-11-15 19:44:54 UTC (rev 5270)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_7.f 2006-11-15 20:04:40 UTC (rev 5271)
@@ -30,14 +30,26 @@
c... Program segment to define material model 7:
c
c Model number: 7
-c Model name: ??
-c Number material properties: ??
-c Number state variables: ??
-c Tangent matrix varies with state: ??
-c Material properties: ??
-c ??
-c ??
+c Model name: IsotropicLinearGenMaxwellViscoelastic
+c Number material properties: 9
+c Number state variables: 15
+c Tangent matrix varies with state: True
+c Material properties: Density
+c Young's Modulus
+c Poisson's Ratio
+c Shear Ratio 1
+c Viscosity 1
+c Shear Ratio 2
+c Viscosity 2
+c Shear Ratio 3
+c Viscosity 3
c
+c... Usage notes:
+c For a simple Maxwell model, set one shear ratio to 1.0 and the
+c other two to 0.0.
+c For a standard linear solid, set two shear ratios to 0.0 and the
+c final one to whatever fraction is appropriate.
+c
subroutine mat_prt_7(prop,nprop,matnum,idout,idsk,kw,kp,
& ierr,errstrng)
c
@@ -61,20 +73,68 @@
c
c... local constants
c
- character labelp(3)*15,modelname*24
+ character labelp(9)*15,modelname*49
data labelp/"Density",
& "Young's modulus",
- & "Poisson's ratio"/
- parameter(modelname="???????")
+ & "Poisson's ratio",
+ & "Shear Ratio 1",
+ & "Viscosity 1",
+ & "Shear Ratio 2",
+ & "Viscosity 2",
+ & "Shear Ratio 3",
+ & "Viscosity 3"/
integer mattype
parameter(mattype=7)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="mat_prt_7"
+ integer i
+ double precision sfrac
c
+cdebug write(6,*) "Hello from mat_prt_5_f!"
+c
+ ierr=0
+ modelname="Isotropic Linear Generalized Maxwell Viscoelastic"
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="mat_prt_5"
+ return
+ end if
+c
+c... output plot results
+c
+ if(idsk.eq.ione) then
+ write(kp,"(3i7)",err=10) matnum,mattype,nprop
+ write(kp,"(1pe15.8,20(2x,1pe15.8))",err=10) (prop(i),i=1,nprop)
+ else if(idsk.eq.itwo) then
+ write(kp,err=10) matnum,mattype,nprop
+ write(kp,err=10) prop
+ end if
+c
+c... output ascii results, if desired
+c
+ if(idout.gt.izero) then
+ write(kw,700,err=10) matnum,modelname,nprop
+ do i=1,nprop
+ write(kw,710,err=10) labelp(i),prop(i)
+ end do
+ end if
+c
return
+c
+c... error writing to output file
+c
+ 10 continue
+ ierr=4
+ errstrng="mat_prt_7"
+ return
+c
+ 700 format(/,5x,"Material number: ",i7,/,5x,
+ & "Material type: ",a37,/,5x,
+ & "Number of properties: ",i7,/)
+ 710 format(15x,a15,3x,1pe15.8)
+c
end
c
c
@@ -100,10 +160,30 @@
character errstrng*(*)
double precision dmat(nddmat),prop(nprop)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="elas_mat_7"
+ integer i,j
+ double precision e,pr,pr1,pr2,pr3,fac,dd,od,ss
+c
+cdebug write(6,*) "Hello from elas_mat_7_f!"
+c
+ call fill(dmat,zero,nddmat)
+ e=prop(2)
+ pr=prop(3)
+ pr1=one-pr
+ pr2=one+pr
+ pr3=one-two*pr
+ fac=e/(pr2*pr3)
+ dd=pr1*fac
+ od=pr*fac
+ ss=half*pr3*fac
+ do i=1,3
+ dmat(iddmat(i,i))=dd
+ dmat(iddmat(i+3,i+3))=ss
+ do j=i+1,3
+ dmat(iddmat(i,j))=od
+ end do
+ end do
return
end
c
@@ -111,10 +191,21 @@
subroutine elas_strs_7(prop,nprop,state,state0,ee,scur,dmat,tmax,
& nstate,nstate0,ierr,errstrng)
c
-c... subroutine to compute stresses for the elastic solution.
+c... subroutine to compute stresses for the elastic solution. For this
+c material, there are 3 sets of state variables: total stress,
+c total strain, and a viscous state variable for each Maxwell
+c element. The minimum Maxwell time is computed, even though this
+c is the elastic solution, as an aid in determining the proper time
+c step size for the next step.
c The current total strain is contained in ee and the computed
-c total strerss should be copied to scur.
+c total stress should be copied to scur.
c
+c state(1:6) = Cauchy stress
+c state(7:12) = linear strain
+c state(13:15) = viscous state variable
+c
+c The state0 array contains initial stresses.
+c
include "implicit.inc"
c
c... parameter definitions
@@ -131,14 +222,92 @@
double precision dmat(nddmat),tmax
character errstrng*(*)
c
-c... return error code, as this material is not yet defined
+c... local variables
c
- ierr=101
- errstrng="elas_strs_7"
+ integer i
+ double precision sfrac,e,pr,vise,rmu,rmue,tmaxe
+c
+cdebug write(6,*) "Hello from elas_strs_7_f!"
+c
+ ierr=0
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="mat_prt_5"
+ return
+ end if
+c
+c... compute elastic stresses assuming material matrix has already
+c been computed.
+c
+ call dcopy(nstr,ee,ione,state(7),ione)
+ call dcopy(nstr,state0,ione,state,ione)
+ call dspmv("u",nstr,one,dmat,state(7),ione,one,state,ione)
+ call dcopy(nstr,state,ione,scur,ione)
+c
+c... compute Maxwell time for current stress state
+c
+ e=prop(2)
+ pr=prop(3)
+ rmu=half*e/(one+pr)
+ tmax=big
+ do i=1,3
+ rmue=prop(2*(i-1)+4)*rmu
+ vise=prop(2*(i-1)+5)
+ if(rmue.ne.zero) tmax=min(tmax,vise/rmue)
+ end do
return
end
c
c
+ function compdq(deltp,vis,rmu)
+c
+c... function to compute the viscous state variable increment delta-q.
+c
+ include "implicit.inc"
+c
+c... parameter definitions
+c
+ include "ndimens.inc"
+ include "rconsts.inc"
+ integer nterms
+ double precision tfrac
+ parameter(nterms=5,tfrac=1.0d-5)
+c
+c... subroutine arguments
+c
+ double precision compdq,deltp,vis,rmu
+c
+c... intrinsic functions
+c
+ intrinsic exp
+c
+c... local variables
+c
+ integer i
+ double precision tau,fsign,fac,frac
+c
+ compdq=zero
+ if(rmu.eq.zero) return
+ tau=vis/rmu
+ if(tau.lt.tfrac*deltp) then
+ fsign=one
+ fac=one
+ frac=one
+ compdq=one
+ do i=2,nterms
+ fac=fac*dble(i)
+ fsign=-one*fsign
+ frac=frac*(deltp/tau)
+ compdq=compdq+fsign*frac/fac
+ end do
+ else
+ compdq=tau*(one-exp(-deltp/tau))/deltp
+ end if
+ return
+ end
+c
+c
subroutine td_matinit_7(state,dstate,state0,dmat,prop,rtimdat,
& rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
& errstrng)
@@ -174,16 +343,57 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... external functions
+c
+ double precision compdq
+ external compdq
+c
+c... local variables
+c
+ double precision sfrac,frac0,e,pr,rmu,rmue,vise,bulkm,sfac,shfac
+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_matinit_7_f!"
c
- ierr=101
- errstrng="td_matinit_7"
+ ierr=0
+ sfrac=prop(4)+prop(6)+prop(8)
+ if(sfrac.lt.0.0d0.or.sfrac.gt.1.0d0) then
+ ierr=116
+ errstrng="mat_prt_5"
+ return
+ end if
+c
+ call fill(dmat,zero,nddmat)
+ e=prop(2)
+ pr=prop(3)
+ rmu=half*e/(one+pr)
+ bulkm=third*e/(one-two*pr)
+ frac0=one-sfrac
+ shfac=frac0
+ tmax=big
+ do i=1,3
+ rmue=prop(2*(i-1)+4)
+ vise=prop(2*(i-1)+5)
+ if(rmue.ne.zero) then
+ tmax=min(tmax,vise/rmue)
+ shfac=shfac+rmue*compdq(deltp,vise,rmue)
+ end if
+ end do
+ shfac=third*rmu*shfac
+ dmat(iddmat(1,1))=bulkm+two*shfac
+ dmat(iddmat(2,2))=dmat(iddmat(1,1))
+ dmat(iddmat(3,3))=dmat(iddmat(1,1))
+ dmat(iddmat(1,2))=bulkm-shfac
+ dmat(iddmat(1,3))=dmat(iddmat(1,2))
+ dmat(iddmat(2,3))=dmat(iddmat(1,2))
+ dmat(iddmat(4,4))=half*three*shfac
+ dmat(iddmat(5,5))=dmat(iddmat(4,4))
+ dmat(iddmat(6,6))=dmat(iddmat(4,4))
return
end
c
More information about the cig-commits
mailing list