Skip to content
Snippets Groups Projects
Commit ea656f23 authored by Dave Whipp's avatar Dave Whipp
Browse files

Added power-law viscosity scaling factor

parent 640903c9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/make_cut.f90
+ 11
− 5
View file @ ea656f23
......@@ -251,10 +251,10 @@ if (icut.eq.0) then
mat(matel)%plasticity_type, &
mat(matel)%plasticity_ss_type_coh, &
mat(matel)%plasticity_ss_type_phi, &
mat(matel)%plasticity_parameters,u,v,w,temp,pressure,strain,&
e2dp,is_plastic_temp,nnode,f,r0,s0,t0,rst,ileaves, &
eviscosityp,vbounded_temp,yield_ratio_temp,frict_angle_temp,&
threadinfo,weight)
mat(matel)%plasticity_parameters,mat(matel)%fviscosity,u,v, &
w,temp,pressure,strain,e2dp,is_plastic_temp,nnode,f,r0,s0, &
t0,rst,ileaves,eviscosityp,vbounded_temp,yield_ratio_temp, &
frict_angle_temp,threadinfo,weight)
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
......@@ -372,6 +372,7 @@ double precision :: eps,viscosity,density,penal,expon,activ,expan,diffusivity
double precision :: heat,eviscosityp,weight,yield_ratio_temp,frict_angle_temp
double precision :: plasticity_parameters(9)
double precision :: aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision :: fvisc
integer :: matel,i
character(len=8) :: plasticity_type
character(len=16) :: plasticity_ss_type_coh,plasticity_ss_type_phi
......@@ -404,6 +405,7 @@ case (1)
plasticity_ss_type_coh=mat(matel)%plasticity_ss_type_coh
plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
plasticity_parameters=mat(matel)%plasticity_parameters
fvisc=mat(matel)%fviscosity
case(2) ! Assign material properties of volumetric minority material
matel=params%materialn(sum(minloc(vol_lsf)))
......@@ -424,6 +426,7 @@ case(2)
plasticity_ss_type_coh=mat(matel)%plasticity_ss_type_coh
plasticity_ss_type_phi=mat(matel)%plasticity_ss_type_phi
plasticity_parameters=mat(matel)%plasticity_parameters
fvisc=mat(matel)%fviscosity
case default
if (matrule.ne.0) write(*,*) 'Invalid matrule value, using divFEM'
......@@ -435,6 +438,7 @@ case default
expan=0.d0
diffusivity=0.d0
heat=0.d0
fvisc=0.d0
do i=1,nlsf
matel=params%materialn(i)
! Check for overriding material transitions
......@@ -450,6 +454,7 @@ case default
expan=expan+vol_lsf(i)*mat(matel)%expansion
diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
heat=heat+vol_lsf(i)*mat(matel)%heat
fvisc=fvisc+vol_lsf(i)*mat(matel)%fviscosity
enddo
matel=params%materialn(0)
if (vol_lsf0.ge.0.015625-eps) matnum=matnum+100
......@@ -460,6 +465,7 @@ case default
activ=activ+vol_lsf0*mat(matel)%activationenergy
expan=expan+vol_lsf0*mat(matel)%expansion
heat=heat+vol_lsf0*mat(matel)%heat
fvisc=fvisc+vol_lsf0*mat(matel)%fviscosity
diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
diffusivity=1.d0/diffusivity ! Note that some properties add geometrically not algebraically
matel=params%materialn(sum(maxloc(vol_lsf)))
......@@ -481,7 +487,7 @@ endif
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,viscosity,density, &
penal,expon,activ,expan,diffusivity,heat,plasticity_type, &
plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,u,v,w,temp,pressure,strain,e2dp, &
plasticity_parameters,fvisc,u,v,w,temp,pressure,strain,e2dp, &
is_plastic_temp,nnode,f,r0,s0,t0,rst,ileaves,eviscosityp, &
vbounded_temp,yield_ratio_temp,frict_angle_temp,threadinfo, &
weight)
......
src/make_matrix.f90
+ 10
− 7
View file @ ea656f23
......@@ -20,10 +20,10 @@ subroutine make_matrix (params,ndof,ael,bel,icon,xg,yg,zg,kfix,viscosity0, &
density,penalty,expon,activationenergy,expansion, &
diffusivity,heat,plasticity_type, &
plasticity_ss_type_coh,plasticity_ss_type_phi, &
plasticity_parameters,unode,vnode,wnode,temp,pressure, &
strain,e2dp,is_plastic,nnode,ffff,r0,s0,t0,rst,ileaves,&
eviscosity,vbounded,yield_ratio,frict_angle,threadinfo,&
weight)
plasticity_parameters,fviscosity,unode,vnode,wnode, &
temp,pressure,strain,e2dp,is_plastic,nnode,ffff,r0,s0, &
t0,rst,ileaves,eviscosity,vbounded,yield_ratio, &
frict_angle,threadinfo,weight)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
......@@ -85,6 +85,7 @@ double precision heat
character (len=8) plasticity_type
character (len=16) plasticity_ss_type_coh,plasticity_ss_type_phi
double precision plasticity_parameters(9)
double precision :: fviscosity
double precision unode(nnode),vnode(nnode),wnode(nnode)
double precision temp(nnode)
double precision pressure
......@@ -236,7 +237,7 @@ do iint=1,nint
J3d=third_invariant (exx,eyy,ezz,exy,eyz,ezx)
e2d=0.d0
if (J2d.ne.0.d0) e2d=sqrt(J2d)
if (J2d.gt.eps) e2d=sqrt(J2d)
straintot=0.d0
do k=1,mpe
......@@ -253,8 +254,10 @@ do iint=1,nint
e2dref=1.d0
viscosity=viscosity0
if (e2d.ne.0.d0 .and. expon.ne.1.d0) viscosity=viscosity*(e2d/e2dref)**(1.d0/expon-1.d0)
viscosity=viscosity*exp(activationenergy/(expon*Rgas*(temperature*params%tempscale+273.15d0)))
if (e2d.gt.eps) then
if (expon.gt.1.d0+eps .or. expon.lt.1.d0-eps) viscosity=viscosity*(e2d/e2dref)**(1.d0/expon-1.d0)
endif
viscosity=viscosity*fviscosity*exp(activationenergy/(expon*Rgas*(temperature*params%tempscale+273.15d0)))
friction_angle=0.d0
if (trim(plasticity_type)/='No') &
......
src/module_definitions.f90
+ 3
− 1
View file @ ea656f23
......@@ -115,13 +115,14 @@ module definitions
! plasticity_parameters are the nine plasticity parameters
! mattrans are the xmin,xmax,ymin,ymax,zmin and zmax for material transitions
! transnum are the corresponding material numbers for the transition
! fviscosity is a scaling factor for power-law viscous behavior
type material
character(len=8) plasticity_type
character(len=16) plasticity_ss_type_coh,plasticity_ss_type_phi
double precision :: plasticity_parameters(9),mattrans(6)
double precision :: density,viscosity,penalty,expon,activationenergy
double precision :: diffusivity,heat,expansion
double precision :: diffusivity,heat,expansion,fviscosity
integer :: transnum(6)
end type material
......@@ -268,6 +269,7 @@ module definitions
integer verminor
integer verstat
integer verrev
integer svnrev
integer dommat
integer,dimension(:),pointer::materialn
double precision ztemp
......
src/read_input_file.f90
+ 6
− 0
View file @ ea656f23
......@@ -244,6 +244,10 @@ do i=0,params%nmat
if (iproc==0) call scanfile (params%infile,'heat'//cm(il:3),mat(i)%heat,ires)
call mpi_bcast(mat(i)%heat,1,mpi_double_precision,0,mpi_comm_world,ierr)
mat(i)%fviscosity=1.d0
if (iproc==0) call scanfile (params%infile,'fviscosity'//cm(il:3),mat(i)%fviscosity,ires)
call mpi_bcast(mat(i)%fviscosity,1,mpi_double_precision,0,mpi_comm_world,ierr)
mat(i)%plasticity_type='No'
if (iproc==0) call scanfile (params%infile,'plasticity_type'//cm(il:3),mat(i)%plasticity_type,ires)
call mpi_bcast(mat(i)%plasticity_type,8,mpi_character,0,mpi_comm_world,ierr)
......@@ -1309,6 +1313,7 @@ if (params%debug.gt.0 .and. iproc.eq.0) then
write(*,'(a,e11.4)') shift//'viscosity ',mat(i)%viscosity
write(*,'(a,e11.4)') shift//'penalty ',mat(i)%penalty
write(*,'(a,e11.4)') shift//'expon ',mat(i)%expon
write(*,'(a,e11.4)') shift//'fviscosity ',mat(i)%fviscosity
write(*,'(a,e11.4)') shift//'diffusivity ',mat(i)%diffusivity
write(*,'(a,e11.4)') shift//'heat ',mat(i)%heat
write(*,'(a,e11.4)') shift//'activationenergy ',mat(i)%activationenergy
......@@ -1579,6 +1584,7 @@ if (params%debug.gt.1) then
write(threadinfo%Logunit,'(a32,e11.4)') 'viscosity ',mat(i)%viscosity
write(threadinfo%Logunit,'(a32,e11.4)') 'penalty ',mat(i)%penalty
write(threadinfo%Logunit,'(a32,e11.4)') 'expon ',mat(i)%expon
write(threadinfo%Logunit,'(a32,e11.4)') 'fviscosity ',mat(i)%fviscosity
write(threadinfo%Logunit,'(a32,e11.4)') 'diffusivity ',mat(i)%diffusivity
write(threadinfo%Logunit,'(a32,e11.4)') 'heat ',mat(i)%heat
write(threadinfo%Logunit,'(a32,e11.4)') 'activationenergy ',mat(i)%activationenergy
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment