[cig-commits] r5367 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Nov 28 14:14:34 PST 2006
Author: willic3
Date: 2006-11-28 14:14:34 -0800 (Tue, 28 Nov 2006)
New Revision: 5367
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f
Log:
Initial version of linear Maxwell model that uses Prony series approach
rather than ESF. This hasn't been tested yet and I may revert back to
the previous version, depending on how tests go.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f 2006-11-28 22:00:16 UTC (rev 5366)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/mat_5.f 2006-11-28 22:14:34 UTC (rev 5367)
@@ -32,13 +32,17 @@
c Model number: 5
c Model name: IsotropicLinearMaxwellViscoelastic
c Number material properties: 4
-c Number state variables: 3
+c Number state variables: 18
c Tangent matrix varies with state: True
c Material properties: Density
-c Young's modulus
-c Poisson's ratio
+c Young's Modulus
+c Poisson's Ratio
c Viscosity
c
+c... This is a revised version of a linear Maxwell model based on the
+c Zienkiewicz and Taylor generalized Maxwell formulation rather than
+c the Bathe effective stress function approach.
+c
subroutine mat_prt_5(prop,nprop,matnum,idout,idsk,kw,kp,
& ierr,errstrng)
c
@@ -46,6 +50,7 @@
c
c Error codes:
c 4: Write error
+c 101: Attempt to use undefined material model
c
include "implicit.inc"
c
@@ -61,12 +66,11 @@
c
c... local constants
c
- character labelp(4)*15,modelname*37
+ character labelp(4)*15,modelname*49
data labelp/"Density",
& "Young's modulus",
& "Poisson's ratio",
& "Viscosity"/
- parameter(modelname="Isotropic Linear Maxwell Viscoelastic")
integer mattype
parameter(mattype=5)
c
@@ -76,6 +80,9 @@
c
cdebug write(6,*) "Hello from mat_prt_5_f!"
c
+ ierr=0
+ modelname="Isotropic Linear Maxwell Viscoelastic"
+c
c... output plot results
c
if(idsk.eq.ione) then
@@ -167,15 +174,17 @@
c
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 viscous strain. The Maxwell time is computed,
-c even though this is the elastic solution, as an aid in determining
-c the proper time step size for the next step.
+c total strain, and a viscous state variable for the 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 stress should be copied to scur.
c
c state(1:6) = Cauchy stress
c state(7:12) = linear strain
-c state(13:18) = viscous strain
+c state(13:18) = dimensionless viscous deviatoric strains for
+c Maxwell element
c
c The state0 array contains initial stresses.
c
@@ -195,28 +204,97 @@
double precision dmat(nddmat),tmax
character errstrng*(*)
c
+c... local constants
+c
+ double precision diag(6),eng(6)
+ data diag/one,one,one,zero,zero,zero/
+ data eng/one,one,one,half,half,half/
+c
c... local variables
c
- double precision e,pr,vis,rmu
+ integer i
+ double precision e,pr,vis,rmu,emean
+ double precision edev(6)
c
cdebug write(6,*) "Hello from elas_strs_5_f!"
c
+ ierr=0
+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 strain quantities to initialize viscous strain
+c
+ emean=third*(ee(1)+ee(2)+ee(3))
+ do i=1,nstr
+ edev(i)=ee(i)*eng(i)-emean*diag(i)
+ end do
+ call dcopy(nstr,edev,ione,state(13),ione)
+c
c... compute Maxwell time for current stress state
c
e=prop(2)
pr=prop(3)
+ rmu=half*e/(one+pr)
vis=prop(4)
- rmu=half*e/(one+pr)
- tmax=vis/rmu
+ tmax=big
+ if(rmu.ne.zero) tmax=min(tmax,vis/rmu)
return
end
c
c
+ function compdq_5(deltp,tau)
+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_5,deltp,tau
+c
+c... intrinsic functions
+c
+ intrinsic exp
+c
+c... local variables
+c
+ integer i
+ double precision fsign,fac,frac
+c
+ compdq_5=zero
+ if(tau.eq.zero) return
+ if(tau.lt.tfrac*deltp) then
+ fsign=one
+ fac=one
+ frac=one
+ compdq_5=one
+ do i=2,nterms
+ fac=fac*dble(i)
+ fsign=-one*fsign
+ frac=frac*(deltp/tau)
+ compdq_5=compdq_5+fsign*frac/fac
+ end do
+ else
+ compdq_5=tau*(one-exp(-deltp/tau))/deltp
+ end if
+ return
+ end
+c
+c
subroutine td_matinit_5(state,dstate,state0,dmat,prop,rtimdat,
& rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,
& errstrng)
@@ -227,9 +305,6 @@
c computed. Thus, for some time-dependent materials, the material
c matrix will only be an approximation of the material matrix for
c the current iteration.
-c As this is a linear viscoelastic material, the tangent matrix is
-c actually independent of the current stresses and strains, so they
-c are not used.
c Note that only the upper triangle is used (or available), as
c dmat is assumed to always be symmetric.
c
@@ -255,9 +330,16 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... external functions
+c
+ double precision compdq_5
+ external compdq_5
+c
c... local variables
c
- double precision e,pr,vis,f1,f2,rmu
+ integer i
+ double precision e,pr,rmu,vis,bulkm,sfac,shfac
+ double precision rmut,tau
c
c... included variable definitions
c
@@ -267,21 +349,38 @@
c
cdebug write(6,*) "Hello from td_matinit_5_f!"
c
+ ierr=0
+c
+c... get material properties
+c
call fill(dmat,zero,nddmat)
e=prop(2)
pr=prop(3)
+ rmu=half*e/(one+pr)
vis=prop(4)
- rmu=half*e/(one+pr)
- tmax=vis/rmu
- f1=third*e/(one-two*pr)
- f2=third/((one+pr)/e+half*alfap*deltp/vis)
- dmat(iddmat(1,1))=f1+two*f2
+ bulkm=third*e/(one-two*pr)
+c
+c... compute Prony series term
+c
+ shfac=one
+ tmax=big
+ tau=zero
+ if(rmu.ne.zero) then
+ tau=vis/rmu
+ tmax=min(tmax,tau)
+ shfac=shfac+compdq_5(deltp,tau)
+ end if
+c
+c... compute stress-strain matrix
+c
+ shfac=third*rmu*shfac
+ dmat(iddmat(1,1))=bulkm+four*shfac
dmat(iddmat(2,2))=dmat(iddmat(1,1))
dmat(iddmat(3,3))=dmat(iddmat(1,1))
- dmat(iddmat(1,2))=f1-f2
+ dmat(iddmat(1,2))=bulkm-two*shfac
dmat(iddmat(1,3))=dmat(iddmat(1,2))
dmat(iddmat(2,3))=dmat(iddmat(1,2))
- dmat(iddmat(4,4))=half*three*f2
+ dmat(iddmat(4,4))=three*shfac
dmat(iddmat(5,5))=dmat(iddmat(4,4))
dmat(iddmat(6,6))=dmat(iddmat(4,4))
return
@@ -292,8 +391,10 @@
& rtimdat,rgiter,ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,
& ierr,errstrng)
c
-c... subroutine to compute the current values for stress, total strain,
-c and viscous strain increment.
+c... subroutine to compute the current stress for the time-dependent
+c solution.
+c Note that only the upper triangle is used (or available), as
+c dmat is assumed to always be symmetric.
c
include "implicit.inc"
c
@@ -318,66 +419,68 @@
include "rgiter_dim.inc"
include "ntimdat_dim.inc"
c
+c... external functions
+c
+ double precision compdq_5
+ external compdq_5
+c
c... local constants
c
- double precision diag(6)
+ double precision diag(6),eng(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
- double precision e,pr,vis,rmu,ae,tf
- double precision emeantdt,smeant,smeantdt,smean0
- double precision gamtau,f1,f2,f3,f4
- double precision epptdt,sdevt,sdev0,sdevtdt,sdevtau
+ integer i,j,ind
+ double precision e,pr,rmu,bulkm,vis
+ double precision etracetdt,emeantdt,smeantdt,sdevtdt,dq
+ double precision edevtdt,edevt,emeant,deltae,rtime
c
-cdebug integer idb
-c
c... included variable definitions
c
include "rtimdat_def.inc"
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-cdebug write(6,*) "Hello from td_strs_5!"
+cdebug write(6,*) "Hello from td_strs_5_f!"
c
ierr=0
c
-c... define material properties
+c... define material properties and dilatational stress
c
e=prop(2)
pr=prop(3)
vis=prop(4)
rmu=half*e/(one+pr)
- ae=(one+pr)/e
- tf=deltp*(one-alfap)
- gamtau=half/vis
- f1=one/(ae+alfap*deltp*gamtau)
- f2=gamtau*tf
- f3=e/(one-two*pr)
- f4=deltp*gamtau
- tmax=vis/rmu
+ bulkm=third*e/(one-two*pr)
+ etracetdt=ee(1)+ee(2)+ee(3)
+ emeantdt=third*etracetdt
+ smeantdt=bulkm*etracetdt
+ emeant=third*(state(7)+state(8)+state(9))
c
-c... define stress and strain quantities
+c... compute Prony series term
c
- emeantdt=third*(ee(1)+ee(2)+ee(3))
- smeantdt=emeantdt*f3
- smeant=third*(state(1)+state(2)+state(3))
- smean0=third*(state0(1)+state0(2)+state0(3))
+ tmax=big
+ rtime=zero
+ dq=zero
+ if(rmu.ne.zero) then
+ rtime=vis/rmu
+ tmax=min(tmax,rtime)
+ dq=compdq_5(deltp,rtime)
+ end if
c
c... compute new stresses and store stress and strain values in dstate
+c
do i=1,nstr
- epptdt=eng(i)*ee(i)-diag(i)*emeantdt-state(i+12)
- sdevt=state(i)-diag(i)*smeant
- sdev0=state0(i)-diag(i)*smean0
- sdevtdt=f1*(epptdt-f2*sdevt+ae*sdev0)
- sdevtau=(one-alfap)*sdevt+alfap*sdevtdt
- dstate(i)=sdevtdt+diag(i)*(smeantdt+smean0)
+ edevtdt=eng(i)*ee(i)-diag(i)*emeantdt
+ edevt=eng(i)*state(i+6)-diag(i)*emeant
+ deltae=edevtdt-edevt
+ dstate(i+12)=state(i+12)*exp(-deltp/rtime)+dq*deltae
+ sdevtdt=two*rmu*dstate(i+12)
+ dstate(i)=diag(i)*smeantdt+sdevtdt+state0(i)
+ dstate(i+6)=ee(i)
scur(i)=dstate(i)
- dstate(i+6)=ee(i)
- dstate(i+12)=f4*sdevtau
end do
c
return
@@ -422,7 +525,7 @@
include "rgiter_def.inc"
include "ntimdat_def.inc"
c
-cdebug write(6,*) "Hello from td_strs_mat_5_f!"
+cdebug write(6,*) "Hello from td_strs_mat_5_f!
c
call td_matinit_5(state,dstate,state0,dmat,prop,rtimdat,rgiter,
& ntimdat,iddmat,tmax,nstate,nstate0,nprop,matchg,ierr,errstrng)
@@ -461,7 +564,7 @@
c
c... local variables
c
- double precision ptmp(10)
+ double precision ptmp(20)
c
cdebug write(6,*) "Hello from prestr_mat_5_f!"
c
@@ -473,8 +576,8 @@
call elas_mat_5(dmat,ptmp,iddmat,nprop,ierr,errstrng)
return
end
-c
-c
+c
+c
subroutine get_state_5(state,dstate,sout,nstate)
c
c... routine to transfer state variables into sout
@@ -525,23 +628,17 @@
double precision sub
data sub/-1.0d0/
c
-cdebug integer idb
+c... local variables
c
-cdebug write(6,*) "Hello from update_state_5_f!"
-c
-cdebug write(6,*) "state-begin:",(state(idb),idb=1,nstate)
-cdebug write(6,*) "dstate-begin:",(dstate(idb),idb=1,nstate)
- call daxpy(nstate-6,sub,state,ione,dstate,ione)
+ call daxpy(nstate,sub,state,ione,dstate,ione)
call daxpy(nstate,one,dstate,ione,state,ione)
-cdebug write(6,*) "state-end:",(state(idb),idb=1,nstate)
-cdebug write(6,*) "dstate-end:",(dstate(idb),idb=1,nstate)
c
return
end
c
-c
+c
c version
-c $Id: mat_5.f,v 1.8 2005/04/01 23:07:43 willic3 Exp $
+c $Id: mat_5.f,v 1.5 2005/04/01 23:10:17 willic3 Exp $
c Generated automatically by Fortran77Mill on Tue May 18 14:18:50 2004
More information about the cig-commits
mailing list