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ðÌKp 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[lTCYie¿¼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@CoÍ·éÔÝ output_dt = TOTAL_TIME / 100 ! t@CoÍ·éÔÝ noutput = 0 ! t@CoÍÌÔ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ðÌKp 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ðÌKp enddo enddo c c==== smoothing call smoothing_CSA ! conservative smoothing c### call smoothing_XSPH ! smoothing in XSPH 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) 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 real8(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, & RHO(MAXP),PMASS(MAXP),CS(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 ! ¹ñfe«W c NUMX1 = 640 NUMY1 = 32 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) c if((j .le. NUMY1/4) .or. (j .gt. NUMY13/4)) then c MATER(ip) = 1 c else c MATER(ip) = 2 c endif 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 c==== material 2: object angle_inc = 90.0D0 / 180.0D0 * PI ! üËpx c RHO2 = 2.785D3 ! Al 2024-T351 CS2 = 5.33D3 ! ¹¬ G2 = 27.6D9 ! ¹ñfe«W c NUMX2 = 32 NUMY2 = 32 xadd = DISA * (NUMX1 / 2 - NUMX2 / 2) + DISA*10 * dcos(angle_inc) yadd = DISA * NUMY1 + DISA * 10 * dsin(angle_inc) c do i=1,NUMX2 do j=1,NUMY2 ip = ip + 1 X0(1,ip) = xadd + DISA * (i - 1) X0(2,ip) = yadd + DISA * (j - 1) enddo enddo
c
NUMP = ip
NP2 = NUMP - NP1
c
return
end
c
c
c
c***********************************************************************
c set parameters for each particle
c***********************************************************************
subroutine set_param
c
implicit real8(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)
c
c==== material 1: plate
c pmass1 = RHO1 * DISA ** 2.0D0
pmass1 = RHO1 * XSIZE * YSIZE * 1.0D0 / NUMX1 / NUMY1
do ip=1,NP1
RHO(ip) = RHO1
PMASS(ip) = pmass1
CS(ip) = CS1
enddo
c
c==== material 2: object
pmass2 = RHO2 * DISA ** 2.0D0
c pmass2 = RHO2 * 4.0D0 / 3.0D0 * PI * (0.5D0 * 1.58D-3) ** 3.0D0
do ip=NP1+1,NUMP
RHO(ip) = RHO2
PMASS(ip) = pmass2
CS(ip) = CS2
enddo
c
return
end
c
c
c
c***********************************************************************
c set constitutive matrix for CFRP lamina
c***********************************************************************
subroutine set_constitutive_matrix
c
implicit real8(a-h,o-z)
c
include 'parameter.f'
c
common/composite/ CMAT(3,3,2),MATER(MAXP)
c
c---- material properties
e11 = 139.0D9
e22 = 9.85D9
g12 = 5.25D9
g23 = 3.8D0
nu12 = 0.25D0
nu23 = 0.38D0
c
c---- 0 degree ply: mater=1
mat = 1
CMAT(1,1,mat) = 1.0D0 / e11
CMAT(2,1,mat) = - nu12 / e11
CMAT(3,1,mat) = 0.0D0
CMAT(1,2,mat) = CMAT(2,1,mat)
CMAT(2,2,mat) = 1.0D0 / e22
CMAT(3,2,mat) = 0.0D0
CMAT(1,3,mat) = 0.0D0
CMAT(2,3,mat) = 0.0D0
CMAT(3,3,mat) = 1.0D0 / g12
c
c---- 90 degree ply: mater=2
mat = 2
CMAT(1,1,mat) = 1.0D0 / e22
CMAT(2,1,mat) = - nu23 / e22
CMAT(3,1,mat) = 0.0D0
CMAT(1,2,mat) = CMAT(2,1,mat)
CMAT(2,2,mat) = 1.0D0 / e22
CMAT(3,2,mat) = 0.0D0
CMAT(1,3,mat) = 0.0D0
CMAT(2,3,mat) = 0.0D0
CMAT(3,3,mat) = 1.0D0 / g23
c
c---- constitutive matrix
do mat=1,2
call inversemat(CMAT(:,:,mat),3)
enddo
c
return
end
c
c
c
c***********************************************************************
c define boundary conditions
c***********************************************************************
subroutine definebc
c
implicit real8(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 real8(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 real8(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 real8(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 ! ÂðCNg
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 real8(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 ! ÂðCNg
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 real8(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 ! lIÈjóÌ»è
do ip=1,NUMP
do l=1,NNEIGH0(ip)
if((RIJ1(l,ip) .gt. 2.0D0H) .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 real8(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,
& RHO(MAXP),PMASS(MAXP),CS(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 real8(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 ! Gruneisenp[
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 real8(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)
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 real8(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 * s2.0D0 + 0.75D0 * s3.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 real8(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 real8(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/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)
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
common/damage/ DMG(MAXP)
c
real8 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
mat = (ip - 1) / NP1 + 1
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 if(mat .ne. 1) then
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
SIGM(1:2,ip) = HENSAS(1:2,ip) - pressure
SIGM(3,ip) = HENSAS(3,ip) ! (1.0D0 - DMG(ip))
c
c else ! mat = 1 (CFRP)
c SIGM(1:3,ip) = SIGM(1:3,ip)
c & + MATMUL(CMAT(:,:,MATER(ip)),epsdot(:)) dt_present
c SIGM(1:2,ip) = SIGM(1:2,ip) * (1.0D0 - DMG(ip))
c endif
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 real8(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,
& RHO(MAXP),PMASS(MAXP),CS(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)
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 real8(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,
& RHO(MAXP),PMASS(MAXP),CS(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
real8 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
real8 SSR(4)/ 0.0D0, ! yXVÅqÉ©©éW ! 0
& -0.29289321881345247559915563789515D0, ! -{2-sqrt(2)}/2
& -1.7071067811865475244008443621048D0, ! -{2+sqrt(2)}/2
& -0.33333333333333333333333333333333D0/ ! 1/3
real8 TRK(4)/ 1.0D0, ! qXVÅkÉ©©éW ! 1
& 0.5857864376269049511983112757903D0, ! 2-sqrt(2)
& 3.4142135623730950488016887242097D0, ! 2+sqrt(2)
& 0.0D0/ ! 0
real8 TTR(4)/ 0.0D0, ! qXVÅqÉ©©éW ! 0
& 0.12132034355964257320253308631455D0, ! {3sqrt(2)-4}/2
& -4.1213203435596425732025330863145D0, ! -{3*sqrt(2)+4}/2
& 0.0D0/ ! 0
c
c==== QENbEM@
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 real8(a-h,o-z)
c
include 'parameter.f'
parameter(EPSMAX=0.5D0, ! jfиÝ
& SGM0MAX=1.2D9) ! jfÍ
c
parameter(FLT=1.45D9,
& FLC=-1.45D9,
& FTT=0.26D9,
& FTC=-1.03D9,
& FLTS=0.456D9)
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)
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
common/damage/ DMG(MAXP)
c
do ip=1,NUMP
mat = (ip - 1) / NP1 + 1
c
c if(mat .ne. 1) then
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
c else
c if(MATER(ip) .eq. 1) then
c fxt = FLT
c fxc = FLC
c fyt = FTT
c fyc = FTC
c fs = FLTS
c else ! mater = 2
c fxt = FTT
c fxc = FTC
c fyt = FTT
c fyc = FTC
c fs = FLTS
c endif
c
c!!! if((SIGM(1,ip) .ge. fxt) .or. (SIGM(1,ip) .le. fxc) .or.
c!!! & (SIGM(2,ip) .ge. fyt) .or. (SIGM(2,ip) .le. fyc) .or.
c!!! & (dabs(SIGM(3,ip)) .ge. fs)) then
c!!! DMG(ip) = 1.0D0
c!!! write(4,1000) 'damage,',X0(1:2,ip)
c!!! endif
c
c hoffman = (-SIGM(1,ip)2.0D0 + SIGM(1,ip)SIGM(2,ip)) /fxc/fxt
c & - SIGM(2,ip) ** 2.0D0 / fyc / fyt
c & + SIGM(3,ip) ** 2.0D0 / fs ** 2.0D0
c & + SIGM(1,ip) * (fxc + fxt) / fxc / fxt
c & + SIGM(2,ip) * (fyc + fyt) / fyc / fyt
cc
c if(hoffman .ge. 1.0D0) then
c DMG(ip) = 1.0D0
c write(4,1000) 'damage,',X0(1:2,ip)
c endif
c endif
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 real8(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, ! JCfp[B (Pa)
& PC=0.015D0, ! JCfp[C
& PN=0.34D0, ! JCfp[n
& PM=1.0D0) ! JCfp[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,
& RHO(MAXP),PMASS(MAXP),CS(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)
common/yield/ SGMY(MAXP),EPSP(MAXP),T(MAXP),JYIELD(MAXP)
c
do ip=1,NUMP
c---- calculate temperature
T(ip) = TREF + E(ip) / CP ! ±q·x
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
else
JYIELD(ip) = 0
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 - tstarPM)
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 real8(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)
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 real8(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,
& RHO(MAXP),PMASS(MAXP),CS(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 real8(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,
& 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
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 real8(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 real8(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 real8(a-h,o-z)
c
dimension X(2,),SIGM(3,)
c
write(7,) 'ZONE I=',NUMP,' F=POINT'
write(8,) 'ZONE I=',NUMP,' F=POINT'
c
do i=1,NUMP
write(7,'(5E12.4)') X(1:2,i),SIGM(1:3,i)
c
sgmx = SIGM(1,i)
sgmy = SIGM(2,i)
tauxy = SIGM(3,i)
c
temp = dsqrt(sgmx2.0D0 - 2.0D0 * sgmx * sgmy + sgmy2.0D0
& + 4.0D0 * tauxy2.0D0)
ps1 = 0.5D0 * (sgmx + sgmy + temp)
ps3 = 0.5D0 * (sgmx + sgmy - temp)
tmx = (ps1 - ps3) * 0.5D0
c
if(tauxy .eq. 0.0D0) then
if(ps1 .eq. sgmx) then
dl = 1.0D0
dm = 0.0D0
else
dm = 1.0D0
dl = 0.0D0
endif
c
else
temp = (sgmx - ps1) ** 2.0D0 / tauxy ** 2.0D0 + 1.0D0
dl = 1.0D0 / dsqrt(temp)
dm = (sgmx - ps1) / tauxy * dl
endif
c
write(8,'(7E12.4)') X(1:2,i),dl,dm,ps1,ps3,tmx
enddo
c
return
end
c
c
c
c********************************************************************
c calculate inverse matrix
c***********************************************************************
subroutine inversemat(a,n)
c
implicit real8(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