[通用算法] 求教首次循环后Pforce（1，：）与Pforce（100，：）不同的原因

program main
use 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