求教首次循环后Pforce（1，：）与Pforce（100，：）不同的原因
program mainuse global
loop=2
lfop=1
call initial
write(*,*)'this is ok1'
call kernel
open(26,file = 'peri_2Dcrack.txt')
write(26,*) 'Title = "Peridynamics Simulation"'
write(26,*) 'Variables = "X", "Y", "D", "C"'
do L = 1,loop
if(mod(L,1)==0 .or. L==1)then
write(*,'(a15,i7,a5,f7.5)') "At present: ",l,"t=",l*dt
end if
call peridynamic
call concentration
do i = 1,totnode
acc(i,1) = (pforce(i,1) + bforce(i,1)) / dens
acc(i,2) = (pforce(i,2) + bforce(i,2)) / dens
vel(i,1) = vel(i,1) + acc(i,1) * dt
vel(i,2) = vel(i,2) + acc(i,2) * dt
disp(i,1) = disp(i,1) + vel(i,1) * dt
disp(i,2) = disp(i,2) + vel(i,2) * dt
end do
!printing results to an output file
if(mod(L,lfop) == 0 ) then
write(26,*) 'zone'
do i = 1, totnode
write(26,'(3(e12.5,3x))') coord(i,1)+disp(i,1), coord(i,2)+disp(i,2), dmg(i,1)
enddo
endif
enddo
close(26)
end program main
module global
integer ndivx, ndivy, totnode, nt, maxfam, nnum, cnode, i, j, loop, nbnd,lfop
real (8) length,width, dx, delta, thick, dens, emod, pratio, area, vol, bc,idist1
real (8) sedload1, sedload2, dt, totime, ctime,fac, radij, dforce1, dforce2
real (8) crlength, scr0, pi, tmpdx, tmpvol, tmpcx, tmpcy, tmpux, tmpuy, dmgpar1, dmgpar2
real (8) scx, scy, scr,pmv,mlmb,kb,tmp,na,cmax,cc1,cc2,cc3,cc4,cc5,cc6
real(8),allocatable::coord(:,:), pforce(:,:), pforceold(:,:), bforce(:,:), stendens(:,:)
real(8),allocatable::fncst(:,:), disp(:,:), vel(:,:), velhalfold(:,:), velhalf(:,:)
real(8),allocatable::acc(:,:), massvec(:,:), enddisp(:,:), endtime(:,:), dmg(:,:)
real(8),allocatable::cc(:), cci(:), theta(:), idist(:,:), nlength(:,:)
real(8),allocatable::PK_I(:,:),PK(:),FF(:,:),F(:),detf(:),sigama(:)
integer(4),allocatable:: numfam(:,:), pointfam(:,:), nodefam(:,:), fail(:,:)
end module global
subroutine initial
use global
ndivx=100;ndivy=100;totnode=ndivx*ndivy;maxfam = 50
allocate(coord(totnode,2), pforce(totnode,2), pforceold(totnode,2), bforce(totnode,2), stendens(totnode,2))
allocate(fncst(totnode,2), disp(totnode,2), vel(totnode,2), velhalfold(totnode,2), velhalf(totnode,2))
allocate(acc(totnode,2), massvec(totnode,2), enddisp(nt,1), endtime(nt,1), dmg(totnode,1))
allocate(cc(totnode), cci(totnode), theta(totnode), idist(totnode,totnode), nlength(totnode,maxfam))
allocate(PK_I(totnode,4),PK(4),FF(totnode,4),F(4),detf(totnode),sigama(totnode))
allocate(numfam(totnode,1), pointfam(totnode,1), nodefam(totnode*maxfam,1), fail(totnode,maxfam))
pi = dacos(-1.0d0);length = 10.0d-7;width = 10.0d-7;dx=1.0d-8;delta=3.015*dx;thick=dx
dens = 8000.0d0;emod = 80.0d9;pratio = 0.22d0;area = dx * dx;vol = area * dx
bc = 9.0d0 * emod / (pi * thick * (delta**3));sedload1 = 9.0d0 / 16.0d0 * emod * 1.0d-6
sedload2 = 9.0d0 / 16.0d0 * emod * 1.0d-6;dt =2.0d-14
pmv=8.5d-6;mlmb=500.0d0;kb=1.38d-23;tmp=300.0d0;na=6.02d23;cmax=1.18d4;crlength = 0.2d-6;radij = dx / 2.0d0;scr0 = 0.04d0
do i = 1, totnode*maxfam
nodefam(i,1) = 0
enddo
do i = 1, nt
enddisp(i,1) = 0.0d0
endtime(i,1) = 0.0d0
enddo
do i = 1,totnode
do j = 1,maxfam
fail(i,j) = 1
enddo
enddo
do i = 1,ndivy
do j = 1,ndivx
nnum = nnum + 1
coord(nnum,1)=(-0.5d0*length)+(0.5d0*dx)+(j-1)*dx
coord(nnum,2)=(-0.5d0*width)+(0.5d0*dx)+(i-1)*dx
end do
end do
write(*,*)coord(totnode,1),coord(totnode,2),nnum
end
subroutine kernel
use global
do i = 1,totnode
if (i.eq.1) then
pointfam(i,1) = 1
else
pointfam(i,1) = pointfam(i-1,1) + numfam(i-1,1)
endif
do j = 1,totnode
idist(i,j)=sqrt((coord(j,1) - coord(i,1))**2 + (coord(j,2) - coord(i,2))**2)
if (i.ne.j) then
if(idist(i,j).le.delta) then
numfam(i,1) = numfam(i,1) + 1
nodefam(pointfam(i,1)+numfam(i,1)-1,1) = j
endif
endif
enddo
write(*,*)i
enddo
do i = 1,totnode
do j = 1,numfam(i,1)
cnode = nodefam(pointfam(i,1)+j-1,1)
if ((coord(cnode,2) > 0.0d0).and.(coord(i,2) < 0.0d0)) then
if ((dabs(coord(i,1)) - (crlength / 2.0d0)).le.1.0d-10) then
fail(i,j) = 0
elseif ((dabs(coord(cnode,1)) - (crlength / 2.0d0)).le.1.0d-10) then
fail(i,j) = 0
endif
elseif ((coord(i,2) > 0.0d0).and.(coord(cnode,2) < 0.0d0)) then
if((dabs(coord(i,1)) - (crlength / 2.0d0)).le.1.0d-10) then
fail(i,j) = 0
elseif((dabs(coord(cnode,1)) - (crlength / 2.0d0)).le.1.0e-10) then
fail(i,j) = 0
endif
endif
enddo
enddo
!Loading 1
do i = 1,totnode
disp(i,1) = 0.001d0 * coord(i,1)
disp(i,2) = 0.0d0
enddo
do i = 1,totnode
stendens(i,1) = 0.0d0
do j = 1,numfam(i,1)
cnode = nodefam(pointfam(i,1)+j-1,1)
idist(i,j) = dsqrt((coord(cnode,1) - coord(i,1))**2 + (coord(cnode,2) - coord(i,2))**2)
nlength(i,j) = dsqrt((coord(cnode,1) + disp(cnode,1) - coord(i,1) - disp(i,1))**2 &
+ (coord(cnode,2) + disp(cnode,2) - coord(i,2) - disp(i,2))**2)
if (idist(i,j).le.delta-radij) then
fac = 1.0d0
elseif (idist(i,j).le.delta+radij) then
fac = (delta+radij-idist(i,j))/(2.0d0*radij)
else
fac = 0.0d0
endif
stendens(i,1)=stendens(i,1)+0.5d0*0.5d0*bc*((nlength(i,j)-idist(i,j))/idist(i,j))**2*idist(i,j)*vol*fac
enddo
fncst(i,1) = sedload1 / stendens(i,1)
if(i==101)then
stop
end if
enddo
!Loading 2
do i = 1,totnode
disp(i,1) = 0.0d0
disp(i,2) = 0.001d0 * coord(i,2)
enddo
do i = 1,totnode
stendens(i,2) = 0.0d0
do j = 1,numfam(i,1)
cnode = nodefam(pointfam(i,1)+j-1,1)
idist(i,j) = dsqrt((coord(cnode,1) - coord(i,1))**2 + (coord(cnode,2) - coord(i,2))**2)
nlength(i,j) = dsqrt((coord(cnode,1) + disp(cnode,1)-coord(i,1)-disp(i,1))**2 &
+ (coord(cnode,2) + disp(cnode,2) - coord(i,2) - disp(i,2))**2)
if (idist(i,j).le.delta-radij) then
fac = 1.0d0
elseif (idist(i,j).le.delta+radij) then
fac = (delta+radij-idist(i,j))/(2.0d0*radij)
else
fac = 0.0d0
endif
stendens(i,2)=stendens(i,2)+0.5d0*0.5d0*bc*((nlength(i,j)-idist(i,j))/idist(i,j))**2*idist(i,j)*vol*fac
enddo
fncst(i,2) = sedload2 / stendens(i,2)
enddo
stop
do i = 1,totnode
vel(i,1) = 0.0d0
vel(i,2) = 0.0d0
disp(i,1) = 0.0d0
disp(i,2) = 0.0d0
enddo
end
subroutine peridynamic
use global
open(24,file='pforce.txt')
open(25,file='pforce1.txt')
write(24,*) 'Variables = "X", "Y", "pforce1", "pforce2","dmg"'
do i = 1,totnode
dmgpar1 = 0.0d0
dmgpar2 = 0.0d0
pforce(i,1) = 0.0d0
pforce(i,2) = 0.0d0
do j = 1,numfam(i,1)
cnode=nodefam(pointfam(i,1)+j-1,1)
idist(i,j) = dsqrt((coord(cnode,1) - coord(i,1))**2 + (coord(cnode,2) - coord(i,2))**2)
nlength(i,j)=dsqrt((coord(cnode,1)+disp(cnode,1)-(coord(i,1)+disp(i,1)))**2 &
+(coord(cnode,2)+disp(cnode,2)-(coord(i,2)+disp(i,2)))**2)
!Volume correction
if (idist(i,j).le.delta-radij) then
fac = 1.0d0
elseif (idist(i,j).le.delta+radij) then
fac = (delta+radij-idist(i,j))/(2.0d0*radij)
else
fac = 0.0d0
endif
if (dabs(coord(cnode,2) - coord(i,2)).le.1.0d-10) then
theta(i) = 0.0d0
elseif (dabs(coord(cnode,1) - coord(i,1)).le.1.0d-10) then
theta(i) = 90.0d0 * pi / 180.0d0
else
theta(i) = datan(dabs(coord(cnode,2) - coord(i,2)) / dabs(coord(cnode,1) - coord(i,1)))
endif
!Determination of the surface correction between two material points
scx = (fncst(i,1) + fncst(cnode,1)) / 2.0d0
scy = (fncst(i,2) + fncst(cnode,2)) / 2.0d0
scr = 1.0d0 / (((dcos(theta(i)))**2 / (scx)**2) + ((dsin(theta(i)))**2 / (scy)**2))
scr = dsqrt(scr)
if (fail(i,j).eq.1) then
dforce1=bc*((nlength(i,j)-idist(i,j))/idist(i,j)-pmv*(cc(i)+cc(cnode))/2)*vol*scr*fac*&
(coord(cnode,1)+disp(cnode,1)-coord(i,1)-disp(i,1))/nlength(i,j)
dforce2=bc*((nlength(i,j)-idist(i,j))/idist(i,j)-pmv*(cc(i)+cc(cnode))/2)*vol*scr*fac*&
(coord(cnode,2)+disp(cnode,2)-coord(i,2)-disp(i,2))/nlength(i,j)
pforce(i,1) = pforce(i,1) + dforce1
pforce(i,2) = pforce(i,2) + dforce2
if(i==101)then
stop
end if
endif
!Definition of a no-fail zone
if (dabs((nlength(i,j) - idist(i,j)) / idist(i,j)) > scr0) then
if (dabs(coord(i,2)).le.(length/4.0d0)) then
fail(i,j) = 0
endif
endif
dmgpar1 = dmgpar1 + fail(i,j) * vol * fac
dmgpar2 = dmgpar2 + vol * fac
enddo
dmg(i,1) = 1.0d0 - dmgpar1 / dmgpar2
end do
stop
end 调试。v.fcode.cn有调试的视频。
页:
[1]