12 Dec 2013 5:07
#13462
c#######################################################################
c
c SPH: Smoothed Particle Hydrodynamics
c
c#######################################################################
c
program SPH
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
parameter(
& TOTAL_TIME=10.0D-6, ! ƒVƒ~ƒ…ƒŒ[ƒVƒ‡ƒ“‚Ìő厞ŠÔ
& MAXLOOP=100000) ! ƒVƒ~ƒ…ƒŒ[ƒVƒ‡ƒ“‚ÌÅ‘åŒJ‚è•Ô‚µ”
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
c
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
c
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
c
common/composite/ CMAT(3,3,2),MATER(MAXP)
c
common/boundary/
& V_BC_GIVEN(NDIM,MAXBC),V_LOAD_GIVEN(NDIM,MAXBC),
& I_BC_GIVEN(NDIM,MAXBC),I_LOAD_GIVEN(NDIM,MAXBC),
& N_BC_GIVEN(NDIM),N_LOAD_GIVEN(NDIM),SPEED
c
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
c
common/param/ H0,H,DT,PDT
c
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
common/kernel/ GRADW(NDIM,NEIGHBOR,MAXP)
c
common/stress/
& SIGM(3,MAXP),HENSAS(3,MAXP),HENSASOLD(3,MAXP),STRNDOT(3,MAXP),STRN(3,MAXP),STRNOLD(3,MAXP)
c
common/vardot/ DRHODT(MAXP),DVDT(NDIM,MAXP),DEDT(MAXP)
c
common/rkgill/ DX(NDIM,MAXP),DV(NDIM,MAXP),DR(MAXP),DE(MAXP)
c
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
c
common/damage/ DMG(MAXP)
c
common/phase/ MELT(MAXP)
c
c==== output files
open(2,file='rst.txt')
open(3,file='out_damage.plt')
open(4,file='rst_damage.csv')
open(7,file='out_stress.plt')
open(100,file='check_coordinates.csv')
c
c==== input
call inp(angle_inc)
c
call set_param
c
c==== initialize
do ip=1,NUMP
X(1:NDIM,ip) = X0(1:NDIM,ip) ! À•W
E(ip) = 0.0D0 ! ƒGƒlƒ‹ƒM[
do k=1,3
SIGM(k,ip) = 0.0D0 ! ‰ž—Í
HENSAS(k,ip) = 0.0D0 ! •Î·‰ž—Í
HENSASOLD(k,ip) = 0.0D0
STRNDOT(k,ip) = 0.0D0 ! ‚Ђ¸‚Ý‘¬“x
STRN(k,ip) = 0.0D0
STRNOLD(k,ip) = 0.0D0
enddo
c
DX(1:NDIM,ip) = 0.0D0 ! ƒ‹ƒ“ƒQEƒNƒbƒ^–@
DV(1:NDIM,ip) = 0.0D0 ! ‚Å‚ÌŽžŠÔÏ•ª‚ÅŽg‚¤
DR(ip) = 0.0D0
DE(ip) = 0.0D0
c
SGMY(ip) = YIELD0 ! ~•š‰ž—Í
EPSP(ip) = 0.0D0 ! ‘Š“–‘Y«‚Ђ¸‚Ý
T(ip) = 300.0D0 ! ‰·“x
c
DMG(ip) = 0.0D0 ! ‘¹ƒpƒ‰ƒ[ƒ^
MELT(ip) = 0.0D0
enddo
c
c---- define boundary conditions
call definebc
c
V(1:NDIM,1:NUMP) = 0.0D0 ! ‘¬“x‚̉Šú‰»
c
do ip=NP1+1,NUMP ! ”òãđ̂̉‘¬
V(1,ip) = SPEED * dcos(angle_inc)
V(2,ip) = SPEED * dsin(angle_inc)
enddo
c
do k=1,NDIM
do i=1,N_BC_GIVEN(k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i) ! ‹«ŠEðŒ‚Ì“K—p
enddo
enddo
c
c==== output
c### call fileoutput(NUMP,X)
call fileoutput2(NUMP,X,DMG,T,V)
call fileoutput3(NUMP,X,SIGM)
c### stop
c
c==== initial neighboring table
call set_neigh
H = H0 ! ƒJ[ƒlƒ‹ƒTƒCƒYi‰e‹¿”¼Œaj
c
c==== calculate relative distance between particles
call update_neigh
call rijcal(1)
c
c
c==== setting
DT = 1.0D-8 ! ŽžŠÔ‚Ý•ƒ¢ti‰Šú’lj
TIMEDIV = 1.0D0
c
c output_dt = TOTAL_TIME / 1000 ! ƒtƒ@ƒCƒ‹o—Í‚·‚鎞ŠÔ‚Ý
output_dt = TOTAL_TIME / 100 ! ƒtƒ@ƒCƒ‹o—Í‚·‚鎞ŠÔ‚Ý
noutput = 0 ! ƒtƒ@ƒCƒ‹o—Í‚ÌŠÔŠu‚É—˜—p
c
c
c#### time increment ####
time_sim = 0.0D0
it = 0
c
do while(time_sim .lt. TOTAL_TIME)
it = it + 1
time_sim = time_sim + DT
c
c==== out for monitor
write(2,'(A,I7,4E13.5)') 'step',it,time_sim,DT,H,enrg
if(mod(it,100) .eq. 1) write(*,'(A,I7,3(A,E12.5))')
& 'step',it,' time',time_sim,' dt',DT,' H',H
c
do k=1,NDIM
do i=1,N_BC_GIVEN(k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i) ! ‹«ŠEðŒ‚Ì“K—p
enddo
enddo
c
c==== update neighboring list at every 10 steps
call update_neigh
c
c==== calculate artificial viscosity
call viscal
c
c
c<<<< Runge-Kutta method >>>>
do irk=1,4
c==== calculate pressure
call calc_press1 ! Mie-GruneisenF–§“x•Ï‰»C‰·“xŒø‰Ê‘å
c### call calc_press2 ! MurnaghanF–§“x•Ï‰»C‰·“xŒø‰Ê¬
c
c==== calculate gradient of kernel function
call rijcal(0)
call calc_gradf
c
c==== calculate stress
call calc_stress(irk)
c
c==== calculate time-derivative of variables
call timederiv
c
c==== update quantities by Runge-Kutta integration
call rugkut_gill(irk)
c
do k=1,NDIM
do i=1,N_BC_GIVEN(k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i) ! ‹«ŠEðŒ‚Ì“K—p
enddo
enddo
c
c==== smoothing
call smoothing_CSA ! conservative smoothing
c### call smoothing_XSPH ! smoothing in XSPH
c do k=1,NDIM
do i=1,N_BC_GIVEN (k)
V(k,I_BC_GIVEN(k,i)) = V_BC_GIVEN(k,i)
enddo
enddo
c
enddo
c<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>
c
c==== calculate relative distance and judge damage between particles
write(4,*) 'simulation time,',time_sim
call rijcal(2)
call judge_damage
c
c==== calculate yield stress by Johnson-Cook model
call calc_yield_JC(DT)
c
c==== update size of the kernel function
call kernel_size
c
c==== update time step
PDT = DT
call update_dt
c
c==== output
nout = int(time_sim / output_dt)
if(nout .ne. noutput) then
noutput = nout
c### call fileoutput(NUMP,X)
call fileoutput2(NUMP,X,DMG,T,V)
call fileoutput3(NUMP,X,SIGM)
endif
c
c==== store variables
enrg = 0.0D0
do ip=1,NUMP
HENSASOLD(1:3,ip) = HENSAS(1:3,ip)
STRNOLD(1:6,ip) = STRN(1:6,ip)
enrg = 0.5D0 * PMASS(ip) * dot_product(V(:,ip),V(:,ip))
& + E(ip) * PMASS(ip)
enddo
c
c
c###############
c### if(mod(it,10) .eq. 1) call fileoutput2(NUMP,X,DMG,T,V)
c### if(it .eq. 200) stop
if(it .eq. MAXLOOP) stop
c###############
enddo ! time increment
c
call CPU_TIME(time)
write(*,*) 'CPU time = ', time
c
stop
end
c
c
c
c***********************************************************************
c input
c***********************************************************************
subroutine inp(angle_inc)
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/composite/ CMAT(3,3,2),MATER(MAXP)
common/param/ H0,H,DT,PDT
c
c==== initial settings
DISA = 0.03125D-3 ! ‰Šú‚Ì—±ŽqŠÔŠu
H0 = DISA * 1.1D0 ! ‰Šú‚̉e‹¿”¼ŒaFVŠÚ
c H0 = DISA * 1.5D0 ! ‰Šú‚̉e‹¿”¼Œa
c
c==== material 1: plate
RHO1 = 2.78D3 ! Al 2024-T351
CS1 = 5.33D3 ! ‰¹‘¬
G1 = 27.6D9 ! ‚¹‚ñ’f’e«ŒW”
PK1 = 2.0D0
c
NUMX1 = 150
NUMY1 = 100
c
ip = 0 ! —±Žq”Ô†F‰Šú‰»
do i=1,NUMX1
do j=1,NUMY1
ip = ip + 1
X0(1,ip) = DISA * (i - 1)
X0(2,ip) = DISA * (j - 1)
enddo
enddo
c
NP1 = ip
c
XMAX = X0(1,NP1) + DISA * 0.5D0
XMIN = X0(1,1) - DISA * 0.5D0
XSIZE = XMAX - XMIN
YMAX = X0(2,NP1) + DISA * 0.5D0
YMIN = X0(2,1) - DISA * 0.5D0
YSIZE = YMAX - YMIN
c
c==== material 2: object
angle_inc = 90.0D0 / 180.0D0 * PI ! “üŽËŠp“x
c
RHO2 = 2.785D3 ! Al 2024-T351
CS2 = 5.33D3 ! ‰¹‘¬
G2 = 27.6D9 ! ‚¹‚ñ’f’e«ŒW”
PK2 = 2.0D0
c
NUMR2 = 2
NUMX2 = 5
NUMY2 = 3
spacing = 3.0D0*NUMR2*DISA
c
dist = 10.0D0*DISA
c
xadd = DISA*(NUMX1-1) / 2.0D0 - spacing*(NUMX2-1)/2.0D00
yadd = DISA*(NUMY1-1) + dist + NUMR2*DISA
c
do k=1,NUMY2
do l=1,NUMX2
do i=1,NUMR2
n=2* i* PI
c
do j=1,n
th = 2* PI* j/n
ip= ip + 1
c
X0(1,ip) = xadd + DISA*(i-1)*dcos(th)
& + spacing*(l-1)
X0(2,ip) = yadd + DISA*(i-1)*dsin(th)
& + spacing*(k-1)
if((i .eq. 1) .and. (j.eq.1)) go to 10000
enddo
10000 continue
enddo
end do
end do
NUMP = ip
NP2 = NUMP - NP1
c
return
end
c
c
c
c***********************************************************************
c set parameters for each particle
c***********************************************************************
subroutine set_param
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
c
c==== material 1: plate
c pmass1 = RHO1 * DISA ** 2.0D0
do ip=1,NP1
RHO(ip) = RHO1
PMASS(ip) = pmass1
CS(ip) = CS1
PK(ip) = PK2
enddo
c
return
end
c
c
c
c==== material 2: object
pmass2 = RHO2 * DISA ** 2.0D0
c
do ip=NP1+1,NUMP
RHO(ip) = RHO2
PMASS(ip) = pmass2
CS(ip) = CS2
PK(IP) = PK2
enddo
c
return
end
c
c
c
c***********************************************************************
c define boundary conditions
c***********************************************************************
subroutine definebc
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/boundary/
& V_BC_GIVEN(NDIM,MAXBC),V_LOAD_GIVEN(NDIM,MAXBC),
& I_BC_GIVEN(NDIM,MAXBC),I_LOAD_GIVEN(NDIM,MAXBC),
& N_BC_GIVEN(NDIM),N_LOAD_GIVEN(NDIM),SPEED
c
c==== plate
c---- x-direction
n_bc_given(1) = 0
c n_bc_given(1) = 2
c i_bc_given(1,1) = 1
c i_bc_given(1,2) = (NUMX1 - 1) * NUMY1 + 1
c v_bc_given(1,1:2) = 0.0D0 ! —¼‚Í‚µ‚ðŒÅ’è
c
c---- y-direction
n_bc_given(2) = 0
c n_bc_given(2) = 2
c i_bc_given(2,1) = 1
c i_bc_given(2,2) = (NUMX1 - 1) * NUMY1 + 1
c v_bc_given(2,1:2) = 0.0D0 ! —¼‚Í‚µ‚ðŒÅ’è
c
c==== object
SPEED = -5.0D2 ! initial speed
c
c==== force boundary
n_load_given(1) = 0
n_load_given(2) = 0
c
return
end
c
c
c
c***********************************************************************
c output data
c***********************************************************************
subroutine fileoutput(NUMP,X)
c
implicit real*8(a-h,o-z)
c
dimension X(2,*)
c
write(3,*) 'ZONE I=',NUMP,' F=POINT'
c
do i=1,NUMP
write(3,'(2E25.15)') X(1:2,i)
enddo
c
return
end
c
c
c
c***********************************************************************
c set table for neighboring particle
c***********************************************************************
subroutine set_neigh
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
c---- material1
call set_neigh2(1,NP1)
c
c---- material2
call set_neigh2(NP1+1,NUMP)
c
c---- connection
do ip=1,NUMP
do j=1,NNEIGH0(ip)
JCONNECT0(j,ip) = 1
enddo
enddo
c
c---- copy initial data
do i=1,NUMP
NNEIGH(i) = NNEIGH0(i)
do l=1,NNEIGH0(i)
NEIGH(l,i) = NEIGH0(l,i)
enddo
enddo
c
return
end
c
c
c
c***********************************************************************
c set table for neighboring particle
c***********************************************************************
subroutine set_neigh2(NPS,NPE)
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
NNEIGH0(NPS:NPE) = 0
c
h02 = (2.0D0 * H0) ** 2.0D0
c
do i=NPS,NPE
do j=i+1,NPE
r2 = (X0(1,j)-X0(1,i)) ** 2.0D0 + (X0(2,j)-X0(2,i)) ** 2.0D0
if(r2 .gt. h02) goto 110
c
NNEIGH0(i) = NNEIGH0(i) + 1 ! ŒÂ”‚ðƒCƒ“ƒNƒŠƒƒ“ƒg
NEIGH0(NNEIGH0(i),i) = j
c
NNEIGH0(j) = NNEIGH0(j) + 1
NEIGH0(NNEIGH0(j),j) = i
c
110 continue
enddo
enddo
c
return
end
c
c
c
c***********************************************************************
c update table for neighboring particle
c***********************************************************************
subroutine update_neigh
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
NNEIGH(1:NUMP) = NNEIGH0(1:NUMP)
h2 = (2.0D0 * H) ** 2.0D0
c
c---- add neighboring particle after initial data
do i=1,NUMP
do j=i+1,NUMP
r2 = (X(1,j) - X(1,i)) ** 2.0D0 + (X(2,j) - X(2,i)) ** 2.0D0
if(r2 .gt. h2) goto 110
c
do l=1,NNEIGH0(i)
if(j .eq. NEIGH0(l,i)) goto 110
enddo
c
NNEIGH(i) = NNEIGH(i) + 1 ! ŒÂ”‚ðƒCƒ“ƒNƒŠƒƒ“ƒg
NEIGH(NNEIGH(i),i) = j ! NNEIGH0‚æ‚èŒã‚ðã‘‚«
c
NNEIGH(j) = NNEIGH(j) + 1
NEIGH(NNEIGH(j),j) = i
c
110 continue
enddo
enddo
c
return
end
c
c
c
c***********************************************************************
c calculate relative distance between particles
c***********************************************************************
subroutine rijcal(jcal0)
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
do ip=1,NUMP
do l=1,NNEIGH(ip)
jp = NEIGH(l,ip)
r2 = 0.0D0
do k=1,NDIM
xij = X(k,ip) - X(k,jp) ! ip - jp
RIJ(k,l,ip) = xij
r2 = r2 + xij ** 2.0D0
enddo
RIJ2(l,ip) = r2
RIJ1(l,ip) = dsqrt(r2)
enddo
enddo
c
c==== calculate RIJ0 before time increment
if(jcal0 .eq. 1) then
do ip=1,NUMP
do l=1,NNEIGH0(ip)
RIJ0(l,ip) = RIJ1(l,ip)
enddo
enddo
c
elseif(jcal0 .eq. 2) then ! ”’l“I‚È”j‰ó‚Ì”»’è
do ip=1,NUMP
do l=1,NNEIGH0(ip)
if((RIJ1(l,ip) .gt. 2.0D0*H) .and.
& (JCONNECT0(l,ip) .eq. 1)) then
JCONNECT0(l,ip) = 0
jp = NEIGH0(l,ip)
xij = (X0(1,ip) + X0(1,jp)) * 0.5D0
yij = (X0(2,ip) + X0(2,jp)) * 0.5D0
write(4,1000) 'numerical,',xij,yij
endif
enddo
enddo
endif
c
return
1000 format(A,2(E12.5,','))
end
c
c
c
c***********************************************************************
c calculate artificial viscosity
c***********************************************************************
subroutine viscal
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(ZETA=0.1D0,A=2.5D0,B=2.5D0)
c parameter(ZETA=0.1D0,A=1.0D0,B=2.0D0)
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
dimension vij(2)
c
h2 = H ** 2.0D0
c
do ip=1,NUMP
do l=1,NNEIGH(ip)
jp = NEIGH(l,ip)
rhob = (RHO(ip) + RHO(jp)) * 0.5D0
cbar = (CS(ip) + CS(jp)) * 0.5D0
c
vij(1:2) = V(1:2,ip) - V(1:2,jp) ! ip - jp
temp = DOT_PRODUCT(vij,RIJ(:,l,ip)) ! RIJ: ip - jp
c
if(temp .lt. 0) then
pmuij = H * temp / (RIJ2(l,ip) + ZETA * h2)
else
pmuij = 0
endif
c
VISC(l,ip) = -(-A * pmuij * cbar + B * pmuij ** 2.0D0) / rhob
enddo
enddo
c
return
end
c
c
c
c***********************************************************************
c calculate pressure: Mie-Gruneisen equation of state
c large density change with thermal effects
c***********************************************************************
subroutine calc_press1
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,
& RHO(MAXP),PMASS(MAXP),CS(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
c
do ip=1,NUMP
mat = (ip - 1) / NP1 + 1 ! Þ—¿”Ô†
c
c---- Þ—¿‚É‚æ‚éƒpƒ‰ƒ[ƒ^‚ÌÝ’è
if(mat .eq. 1) then
gamma = 2.0D0
rho0 = RHO1
s = 1.338D0
else
gamma = 2.0D0 ! Gruneisenƒpƒ‰ƒ[ƒ^
rho0 = RHO2
s = 1.338D0 ! ÕŒ‚”g‘¬“x‚Æ—±Žq‘¬“x‚ÌüŒ`‹ßŽ—
endif
c
a1 = rho0 * CS(ip) ** 2.0D0
a2 = a1 * (1.0D0 + 2.0D0 * (s - 1.0D0))
a3 = a1 * (2.0D0 * (s - 1.0D0) + 3.0D0 * (s - 1.0D0) ** 2.0D0)
c
eta = RHO(ip) / rho0 - 1.0D0
if(eta .gt. 0) then
ph = a1 * eta + a2 * eta ** 2.0D0 + a3 * eta ** 3.0D0
else
ph = a1 * eta
endif
c
P(ip) = (1.0D0 - 0.5D0 * gamma * eta) * ph
& + gamma * RHO(ip) * E(ip)
enddo
c
return
end
c
c
c
c***********************************************************************
c calculate pressure: Murnaghan equation of state
c small density change and thermal effects are negligible
c***********************************************************************
subroutine calc_press2
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
c
do ip=1,NUMP
mat = (ip - 1) / NP1 + 1 ! Þ—¿”Ô†
c
c---- Þ—¿‚É‚æ‚éƒpƒ‰ƒ[ƒ^‚ÌÝ’è
if(mat .eq. 1) then
rho0 = RHO1
else
rho0 = RHO2
endif
c
P(ip) = CS(ip) ** 2.0D0 * (RHO(ip) -rho0)
enddo
c
return
end
c
c
c
c***********************************************************************
c weight (kernel) function
c***********************************************************************
double precision function weightf(r,h)
c
implicit real*8(a-h,o-z)
c
parameter(D=0.45472840883398667362576479638247D0) ! 10 / 7 / PI
c
s = r / h
coef = D / h ** 2.0D0
c
if(s .le. 1.0D0) then
weightf = coef * (1.0D0 - 1.5D0 * s**2.0D0 + 0.75D0 * s**3.0D0)
elseif(s .le. 2.0D0) then
weightf = coef * 0.25D0 * (2.0D0 - s) ** 3.0D0
else
wieghtf = 0.0D0
endif
c
end
c
c
c
c***********************************************************************
c calculate gradient of weight function
c***********************************************************************
subroutine calc_gradf
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(D=0.45472840883398667362576479638247D0) ! 10 / 7 / PI
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
common/kernel/ GRADW(NDIM,NEIGHBOR,MAXP)
common/param/ H0,H,DT,PDT
c
coef = D / H ** 2.0D0
c
do ip=1,NUMP
do l=1,NNEIGH(ip)
s = RIJ1(l,ip) / H
if(s .le. 1.0D0) then
dwds = coef * (-3.0D0 * s + 2.25D0 * s ** 2.0D0)
elseif(s .le. 2.0D0) then
dwds = -0.75D0 * coef * (2.0D0 - s) ** 2.0D0
else
dwds = 0.0D0
endif
c
temp = dwds / H / RIJ1(l,ip)
GRADW(1:2,l,ip) = RIJ(1:2,l,ip) * temp
enddo
enddo
c
return
end
c
c
c
c***********************************************************************
c calculate stress
c***********************************************************************
subroutine calc_stress(IRK)
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/composite/ CMAT(3,3,2),MATER(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
common/kernel/ GRADW(NDIM,NEIGHBOR,MAXP)
common/stress/
& SIGM(3,MAXP),HENSAS(3,MAXP),HENSASOLD(3,MAXP),STRNDOT(3,MAXP)
& STRN(3,MAXP),STRNOLD(3,MAXP)
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
common/damage/ DMG(MAXP)
common/phase/ MELT(MAXP)
c
real*8 TDIV(4)/ ! Runge-Kutta–@‚Å‚ÌŽžŠÔ‚Ý‚ÌŒW”
& 0.5D0,0.5D0,1.0D0,1.0D0/
c
dimension g(2),vji(2),epsdot(3),rotdot(3),hensasdot(3)
c
dt_present = DT * TDIV(IRK)
c
g(1) = G1
g(2) = G2
c
do ip=1,NUMP
c
c---- strain velocity, rotation velocity
do k=1,3
epsdot(k) = 0.0D0 ! initialize
rotdot(k) = 0.0D0
enddo
c
do l=1,NNEIGH0(ip)
jp = NEIGH0(l,ip)
vji(1:2) = V(1:2,jp) - V(1:2,ip) ! jp - ip
volj = PMASS(jp) / RHO(jp)
c
do k=1,2
epsdot(k) = epsdot(k)
& + volj * vji(k) * GRADW(k,l,ip) * JCONNECT0(l,ip)
ccc rotdot(k) = rotdot(k) + 0.0D0
enddo
c
temp1 = vji(1) * GRADW(2,l,ip) * 0.5D0 * volj *JCONNECT0(l,ip)
temp2 = vji(2) * GRADW(1,l,ip) * 0.5D0 * volj *JCONNECT0(l,ip)
epsdot(3) = epsdot(3) + temp1 + temp2
rotdot(3) = rotdot(3) + temp1 - temp2
enddo
c
STRNDOT(1:3,ip) = epsdot(1:3)
c
c
c---- velocity of deviator stress
taikaku = (epsdot(1) + epsdot(2)) / 3.0D0
w12s12x2 = rotdot(3) * HENSAS(3,ip) * 2.0D0
c
hensasdot(1) = 2.0D0 * g(mat) * (epsdot(1) - taikaku) + w12s12x2
hensasdot(2) = 2.0D0 * g(mat) * (epsdot(2) - taikaku) - w12s12x2
c
w12s22 = rotdot(3) * HENSAS(2,ip)
s11w12 = HENSAS(1,ip) * rotdot(3)
hensasdot(3) = 2.0D0 * g(mat) * epsdot(3) + w12s22 - s11w12
c
c---- deviator stress considering yield stress
do k=1,3
cal = hensasdot(k) * dt_present
HENSAS(k,ip) = HENSASOLD(k,ip) + cal
enddo
c
if(JYIELD(ip) .eq. 1) then
s_invariant2 = (HENSAS(1,ip) ** 2.0D0 + HENSAS(2,ip) ** 2.0D0
& +2.0D0 * HENSAS(3,ip) ** 2.0D0) * 0.5D0
ratio = dsqrt(SGMY(ip) ** 2.0D0 / 3.0D0 / s_invariant2)
HENSAS(1:3,ip) = HENSAS(1:3,ip) * ratio
endif
c
c---- stress
if(P(ip) .ge. 0.0D0) then
pressure = P(ip)
else
pressure = P(ip) * (1.0D0 - DMG(ip))
endif
c
if(MELT(IP).eq.1.0d0) then
SIGM(1:2,ip) = -pressure
SIGM(3,ip) = 0.0D0
c
STRN(1:6,ip) = STRNOLD(1:6,ip)
& + epsdot(1:6)* dt_present
c strn_q = -(STRN(l,ip) + STRN(2,ip) + STRN(3,ip))
c c2 = CS(ip)**2.0D0
c temp = 1.0D0 - PK(ip)*strn_q
c P(ip) = RHO(ip)*c2*strn_q/temp/temp
c
c SIGM(1:2,ip) = - P(ip)
c SIGM(3,ip) = 0.0D0
c
else
SIGM(1:2,ip) = HENSAS(1:2,ip)-pressure
SIGM(3,ip) = HENSAS(3,ip) !* (1.0D0 - DMG(ip))
endif
c
enddo ! ip
c
return
end
c
c
c
c***********************************************************************
c calculate time-derivative of variables
c repulsive artificial force (artificial pressure) is introduced
c***********************************************************************
subroutine timederiv
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(EPS=0.3D0,DN=4.0D0) ! repulsive artificial force
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP)
common/kernel/ GRADW(NDIM,NEIGHBOR,MAXP)
common/stress/
& SIGM(3,MAXP),HENSAS(3,MAXP),HENSASOLD(3,MAXP),STRNDOT(3,MAXP)
& STRN(3,MAXP),STRNOLD(3,MAXP)
common/vardot/ DRHODT(MAXP),DVDT(NDIM,MAXP),DEDT(MAXP)
c
dimension vij(2),ari(2),arj(2),arij(2)
c
c\\\\\\ winv = 1.0D0 / weightf(DISA,H)
c
do ip=1,NUMP
calr = 0.0D0
calv1 = 0.0D0
calv2 = 0.0D0
cale = 0.0D0
rho2inv = 1.0D0 / RHO(ip) ** 2.0D0
c
do l=1,NNEIGH(ip)
jp = NEIGH(l,ip)
vij(1:2) = V(1:2,ip) - V(1:2,jp) ! ip - jp
c
if((l .gt. NNEIGH0(ip)) .or. (JCONNECT0(l,ip) .eq. 0)) then
compression = dmin1(0.0D0,-P(ip)) ! ‚‚Ȃª‚Á‚Ä‚¢‚È‚¢ê‡
s1ip = compression
s2ip = compression
s12ip = 0.0D0
compression = dmin1(0.0D0,-P(jp))
s1jp = compression
s2jp = compression
s12jp = 0.0D0
c
c\\\\\\ ari(1:2) = 0.0D0
c\\\\\\ arj(1:2) = 0.0D0
c
else
s1ip = SIGM(1,ip)
s2ip = SIGM(2,ip)
s12ip = SIGM(3,ip)
s1jp = SIGM(1,jp)
s2jp = SIGM(2,jp)
s12jp = SIGM(3,jp)
c
c\\\\\\ do k=1,2 ! Monaghan 2001
c\\\\\\ ari(k) = dmax1(SIGM(k,ip),0.0D0) * (- EPS)
c\\\\\\ arj(k) = dmax1(SIGM(k,jp),0.0D0) * (- EPS)
c enddo
endif
c
c==== repulsive artificial force
c\\\\\\ arij(1:2) = ari(1:2) + arj(1:2)
c\\\\\\ fn = (weightf(RIJ1(l,ip),H) * winv) ** DN
c
c==== calculation
c
cal11 = (s1ip * rho2inv + s1jp / RHO(jp)**2.0D0 + VISC(l,ip))
c\\\\\\ cal11 = (s1ip * rho2inv + s1jp / RHO(jp)**2.0D0 + VISC(l,ip)
c\\\\\\ & +arij(1) * fn)
& * PMASS(jp) * GRADW(1,l,ip)
temp = (s12ip * rho2inv + s12jp / RHO(jp)**2.0D0)
& * PMASS(jp)
cal12 = temp * GRADW(2,l,ip)
cal21 = temp * GRADW(1,l,ip)
cal22 = (s2ip * rho2inv + s2jp / RHO(jp)**2.0D0 + VISC(l,ip))
c\\\\\\ cal22 = (s2ip * rho2inv + s2jp / RHO(jp)**2.0D0 + VISC(l,ip)
c\\\\\\ & +arij(2) * fn)
& * PMASS(jp) * GRADW(2,l,ip)
c
c==== density
calr = calr
& + PMASS(jp) / RHO(jp) * DOT_PRODUCT(vij,GRADW(:,l,ip))
c
c==== velocity
calv1 = calv1 + cal11 + cal12
calv2 = calv2 + cal21 + cal22
c
c==== energy
cale = cale
& - vij(1) * (cal11 + cal12) - vij(2) * (cal21 + cal22)
enddo ! l
c
DRHODT(ip) = RHO(ip) * calr
DVDT(1,ip) = calv1
DVDT(2,ip) = calv2
DEDT(ip) = 0.5D0 * cale
enddo ! ip
c
return
end
c
c
c
c***********************************************************************
c calculate quantities by Runge-Kutta integration
c***********************************************************************
subroutine rugkut_gill(IRK)
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/param/ H0,H,DT,PDT
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP)
common/kernel/ GRADW(NDIM,NEIGHBOR,MAXP)
common/vardot/ DRHODT(MAXP),DVDT(NDIM,MAXP),DEDT(MAXP)
common/rkgill/ DX(NDIM,MAXP),DV(NDIM,MAXP),DR(MAXP),DE(MAXP)
c
real*8 SRK(4)/ 0.5D0, ! yXV‚Åk‚É‚©‚©‚éŒW” ! 1/2
& 0.29289321881345247559915563789515D0, ! {2-sqrt(2)}/2
& 1.7071067811865475244008443621048D0, ! {2+sqrt(2)}/2
& 0.16666666666666666666666666666667D0/ ! 1/6
real*8 SSR(4)/ 0.0D0, ! yXV‚Åq‚É‚©‚©‚éŒW” ! 0
& -0.29289321881345247559915563789515D0, ! -{2-sqrt(2)}/2
& -1.7071067811865475244008443621048D0, ! -{2+sqrt(2)}/2
& -0.33333333333333333333333333333333D0/ ! 1/3
real*8 TRK(4)/ 1.0D0, ! qXV‚Åk‚É‚©‚©‚éŒW” ! 1
& 0.5857864376269049511983112757903D0, ! 2-sqrt(2)
& 3.4142135623730950488016887242097D0, ! 2+sqrt(2)
& 0.0D0/ ! 0
real*8 TTR(4)/ 0.0D0, ! qXV‚Åq‚É‚©‚©‚éŒW” ! 0
& 0.12132034355964257320253308631455D0, ! {3*sqrt(2)-4}/2
& -4.1213203435596425732025330863145D0, ! -{3*sqrt(2)+4}/2
& 0.0D0/ ! 0
c
c==== ƒ‹ƒ“ƒQEƒNƒbƒ^EƒMƒ‹–@
s1 = SRK(IRK) * DT
s2 = SSR(IRK)
t1 = TRK(IRK) * DT
t2 = TTR(IRK)
c
do ip=1,NUMP
do k=1,2
c---- x
X(k,ip) = X(k,ip) + s1 * V(k,ip) + s2 * DX(k,ip)
DX(k,ip) = t1 * V(k,ip) + t2 * DX(k,ip)
c---- v
V(k,ip) = V(k,ip) + s1 * DVDT(k,ip) + s2 * DV(k,ip)
DV(k,ip) = t1 * DVDT(k,ip) + t2 * DV(k,ip)
enddo
c
c---- rho
RHO(ip) = RHO(ip) + s1 * DRHODT(ip) + s2 * DR(ip)
DR(ip) = t1 * DRHODT(ip) + t2 * DR(ip)
c
c---- E
E(ip) = E(ip) + s1 * DEDT(ip) + s2 * DE(ip)
DE(ip) = t1 * DEDT(ip) + t2 * DE(ip)
enddo
c
return
end
c
c
c
c***********************************************************************
c judge damage between particles
c***********************************************************************
subroutine judge_damage
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(EPSMAX=0.5D0, ! ”j’f‚Ђ¸‚Ý
& SGM0MAX=1.2D9) ! ”j’f‰ž—Í
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/composite/ CMAT(3,3,2),MATER(MAXP)
common/particle/ X(NDIM,MAXP),X0(NDIM,MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/stress/
& SIGM(3,MAXP),HENSAS(3,MAXP),HENSASOLD(3,MAXP),STRNDOT(3,MAXP)
& STRN(3,MAXP),STRNOLD(3,MAXP)
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
common/damage/ DMG(MAXP)
c
do ip=1,NUMP
c
if((dabs(P(ip)) .gt. SGM0MAX) .or. (EPSP(ip) .gt. EPSMAX)) then
DMG(ip) = 1.0D0
write(4,1000) 'damage,',X0(1:2,ip)
endif
c
enddo
c
return
1000 format(A,2(E12.5,','))
end
c
c
c
c***********************************************************************
c calculate temperature
c calculate yield stress by Johnson-Cook
c***********************************************************************
subroutine calc_yield_JC(DT)
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
parameter(TREF=300.0D0, ! ŽQƉ·“x(K)
& CP=0.87D3, ! ’舳”ä”M(J/kg/K)
& TMELT=775.0D0, ! —Z“_(K)
& PB=426D6, ! JCƒ‚ƒfƒ‹ƒpƒ‰ƒ[ƒ^B (Pa)
& PC=0.015D0, ! JCƒ‚ƒfƒ‹ƒpƒ‰ƒ[ƒ^C
& PN=0.34D0, ! JCƒ‚ƒfƒ‹ƒpƒ‰ƒ[ƒ^n
& PM=1.0D0) ! JCƒ‚ƒfƒ‹ƒpƒ‰ƒ[ƒ^m
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/stress/
& SIGM(3,MAXP),HENSAS(3,MAXP),HENSASOLD(3,MAXP),STRNDOT(3,MAXP)
& STRN(3,MAXP),STRNOLD(3,MAXP)
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
common/phase/ MELT(MAXP)
c
do ip=1,NUMP
c---- calculate temperature
T(ip) = TREF + E(ip) / CP ! —±Žq‰·“x
c
c---- judge melt
if(T(ip).ge.TMELT) then
MELT(ip) = 1.0D0
else
Melt(ip)=0.0D0
end if
c
c---- judge yield
s_invariant2 = (HENSAS(1,ip) ** 2.0D0 + HENSAS(2,ip) ** 2.0D0
& +2.0D0 * HENSAS(3,ip) ** 2.0D0) * 0.5D0
if(s_invariant2 .ge. SGMY(ip)**2.0D0 / 3.0D0) then ! ~•š”»’è
JYIELD(ip) = 1
endif
c
c---- if yielding
if(JYIELD(ip) .eq. 1) then
c---- calculate equivalent plastic strain
dk2 = ((STRNDOT(1,ip) - STRNDOT(2,ip)) ** 2.0D0
& +STRNDOT(2,ip) ** 2.0D0 ! eps_z_dot
& +STRNDOT(1,ip) ** 2.0D0 ! eps_z_dot
& +1.5D0 * STRNDOT(3,ip) ** 2.0D0) / 6.0D0
c
epspdot = dsqrt(4.0D0 * dk2 / 3.0D0)
EPSP(ip) = EPSP(ip) + epspdot * DT
c
c---- Johnson-Cook model
tstar = (T(ip) - TREF) / (TMELT - TREF) ! –³ŽŸŒ³‰»‰·“x
c!!! SGMY(ip) = (YIELD0 + PB * EPSP(ip)**PN)
c!!! & * (1.0D0 + PC * dlog(epspdot)) ! eps0dot = 1.0
c!!! & * (1.0D0 - tstar**PM)
c!!! if(SGMY(ip) .lt. 0.0D0) SGMY(ip) = YIELD0
SGMY(ip) = YIELD0
endif
enddo
c
return
end
c
c
c
c***********************************************************************
c update size of the kernel function
c***********************************************************************
subroutine kernel_size
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/param/ H0,H,DT,PDT
common/stress/
& SIGM(3,MAXP),HENSAS(3,MAXP),HENSASOLD(3,MAXP),STRNDOT(3,MAXP),
& STRN(3,MAXP),STRNOLD(3,MAXP)
c
avg_strndot = 0.0D0 ! ‚·‚ׂÄÞ.V < 0‚¾‚Á‚½‚çH0‚Ì‚Ü‚Ü
c
do ip=1,NUMP
temp = STRNDOT(1,ip) + STRNDOT(2,ip)
if(temp .gt. avg_strndot) avg_strndot = temp
enddo
avg_strndot = avg_strndot / 3.0D0
c
dh = H * avg_strndot * DT
H = H + dh
c
return
end
c
c
c
c***********************************************************************
c update time step
c***********************************************************************
subroutine update_dt
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(OMEGA=0.3D0) !,CMAX=0.2D0)
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/param/ H0,H,DT,PDT
c
dtmin = 2.0D-8 ! ‚±‚ê‚æ‚èdt‚ð‘å‚«‚‚µ‚È‚¢
c
do ip=1,NUMP
vp = dsqrt(V(1,ip) ** 2.0D0 + V(2,ip) ** 2.0D0)
temp = OMEGA * H / (CS(ip) + vp) ! CFLðŒ
if(temp .lt. dtmin) dtmin = temp
c
c### temp = DISA * CMAX / vp ! ƒN[ƒ‰ƒ“”—MPS
c### if(temp .lt. dtmin) dtmin = temp
enddo
DT = dtmin
c
return
end
c
c
c
c***********************************************************************
c conservative smoothing: Randles 1996
c***********************************************************************
subroutine smoothing_CSA
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(ACS=0.05D0)
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,PK1,PK2
& RHO(MAXP),PMASS(MAXP),CS(MAXP),PK(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
common/param/ H0,H,DT,PDT
c
do ip=1,NUMP
bunbo = 0.0D0
v1 = 0.0D0
v2 = 0.0D0
density = 0.0D0
energy = 0.0D0
do l=1,NNEIGH0(ip)
jp = NEIGH0(l,ip)
temp = PMASS(jp) * weightf(RIJ1(l,ip),H) / RHO(jp)
& * dble(JCONNECT0(l,ip))
bunbo = bunbo + temp
v1 = v1 + V(1,jp) * temp
v2 = v2 + V(2,jp) * temp
density = density + RHO(jp) * temp
energy = energy + E(jp) * temp
enddo
c
if(bunbo .ne. 0.0D0) then
V(1,ip) = V(1,ip) + ACS * (v1 / bunbo - V(1,ip))
V(2,ip) = V(2,ip) + ACS * (v2 / bunbo - V(2,ip))
RHO(ip) = RHO(ip) + ACS * (density / bunbo - RHO(ip))
E(ip) = E(ip) + ACS * (energy / bunbo - E(ip))
endif
enddo
c
return
end
c
c
c
c***********************************************************************
c conservative smoothing: Monaghan 2001
c***********************************************************************
subroutine smoothing_XSPH
c
implicit real*8(a-h,o-z)
c
include 'parameter.f'
parameter(ZETA=0.5D0)
c
common/model/ DISA,XMAX,XMIN,XSIZE,YMAX,YMIN,YSIZE,
& NUMX1,NUMY1,NUMX2,NUMY2,NUMP,NP1,NP2
common/property/ RHO1,RHO2,CS1,CS2,G1,G2,
& RHO(MAXP),PMASS(MAXP),CS(MAXP)
common/calc/ V(NDIM,MAXP),VISC(NEIGHBOR,MAXP),P(MAXP),E(MAXP)
common/neighbor_table/
& NEIGH(NEIGHBOR,MAXP),NNEIGH(MAXP),
& NEIGH0(NEIGHBOR,MAXP),NNEIGH0(MAXP),JCONNECT0(NEIGHBOR,MAXP),
& RIJ(NDIM,NEIGHBOR,MAXP),RIJ1(NEIGHBOR,MAXP),RIJ2(NEIGHBOR,MAXP),
& RIJ0(NEIGHBOR,MAXP)
common/param/ H0,H,DT,PDT
c
dimension tmp(NDIM)
c
do ip=1,NUMP
tmp(1:2) = 0.0D0
c
do l=1,NNEIGH0(ip)
jp = NEIGH0(l,ip)
rhoavg = (RHO(ip) + RHO(jp)) * 0.5D0
c
do k=1,2
tmp(k) = tmp(k) + PMASS(jp) / rhoavg * (V(k,jp) - V(k,ip))
& * weightf(RIJ1(l,ip),H) * JCONNECT0(l,ip)
enddo
enddo
c
do k=1,2
V(k,ip) = V(k,ip) + ZETA * tmp(k)
enddo
enddo
c
return
end
c
c
c
c***********************************************************************
c output data
c***********************************************************************
subroutine fileoutput2(NUMP,X,VAL1,VAL2,V)
c
implicit real*8(a-h,o-z)
c
dimension X(2,*),VAL1(*),VAL2(*),V(2,*)
c
write(3,*) 'ZONE I=',NUMP,' F=POINT'
c
do i=1,NUMP
write(3,'(6E12.4)') X(1:2,i),VAL1(i),VAL2(i),V(1:2,i)
enddo
c
return
end
c
c
c
c***********************************************************************
c output data
c***********************************************************************
subroutine fileoutput3(NUMP,X,SIGM)
c
implicit real*8(a-h,o-z)
c
dimension X(2,*),SIGM(3,*)
c
write(7,*) 'ZONE I=' ,NUMP, 'F=POINT'
c
do I=1,NUMP
write(7,'(5E12.4)' X(1:2,i),SIGM(1:3,i)
enddo
c
return
end
c
c
c
c***********************************************************************
c calculate inverse matrix
c***********************************************************************
subroutine inversemat(a,n)
c
implicit real*8(a-h,o-z)
c
dimension a(n,*)
c
do k=1,n
aw = a(k,k)
a(k,k) = 1.0d0
do j=1,n
a(k,j) = a(k,j) / aw
enddo
c
do i=1,n
if(i .ne. k) then
aw = a(i,k)
a(i,k) = 0.0d0
do j=1,n
a(i,j) = a(i,j) - aw * a(k,j)
enddo
endif
enddo
enddo
c
return
end
c
c
c
and parameter(include file) [/code] parameter( 1 MAXP=70000, ! ±q 2 NDIM=2, ! ³ 3 NEIGHBOR=100, ! é±qiÉÖW·é±qjÌÅå 4 MAXBC=100) ! «EðÌÅå c parameter(GRAV=9.80665D0) ! m/s^2 c parameter(PI=3.14159265358979323846D0) c parameter(YIELD0=265D6) ! ~Í c
(please help me fixing the error....very urgent thank you)