VDLOAD of LSP

This is a VDLOAD subroutine for Abaqus/Explicit, written using FORTRAN. This code can achieve controllable parameters such as the overlap rate, energy density, and spot shape of the LSP process that can be controlled in actual experiments.

      subroutine vdload (
C Read only (unmodifiable)variables -
     1 nblock, ndim, stepTime, totalTime,
     2 amplitude, curCoords, velocity, dirCos, jltyp, sname,
C Write only (modifiable) variable -
     1 value )
C
      include 'vaba_param.inc'
C
      dimension curCoords(nblock,ndim), velocity(nblock,ndim),
     1  dirCos(nblock,ndim,ndim), value(nblock)
      character*80 sname

C LSP parameters:pulse energy, pulse width, spot size, overlapping ratio
      parameter(pw=9, d=2.d0,ol=5d-1)
      parameter(stpx=3.d0,stpy=2.5d0,xn=7.d0,yn=6.d0,oneperiod=999d-9)
C Shape of pulse laser：timea, timeb, timec, peak pressure
      parameter(timea=0.d0,timeb=9d-9,timec=27d-9,pressure=6.6d3)
C variable in calculation：1，2，3，0.5，pi
      parameter(one=1.d0,two=2.d0,three=3.d0,half=5.d-1,pi=3.1415926d0)
C variable in loops:point (Refers to the serial number of the impact point),ol(overlapping ratio)
      integer ai,aj,ak,point
      double precision al,ys,dp,increase,cx,cy,sx,sy
C Variables about time:point number, start time, endtime, time of one period
      double precision stime,etime,time
C Variables about location:current x, current y, distance,setpoint x,setpoint y,x number,y number,pn
      double precision distance,r
C Variables about energy：uniformed pressure, energy density
      double precision qp,ed
      
      increase=d*(1.d0-ol) 
      dp=floor(totaltime/oneperiod)+one
      point=int(dp)
      r=d/2.0d0
C aj represents the row number，ak refers to directions 1 left 0 right
      aj=((point-1)/xn)+1
      ak=mod(aj,2)
      if(ak==1) then
             al=-1.d0
             sx=stpx
      else
             al=1.d0
             sx=stpx-(xn-1)*increase
      end if
C x coordinate
      sx=sx+al*increase*(point-1-(aj-1)*xn)
c y coordinate
      sy=stpy-increase*(aj-1)
          
      time=totaltime-(point-1)*oneperiod
      do km=1,nblock
          cx=curcoords(km,2)
          cy=-one*curcoords(km,1)
          distance=((cx-sx)**two+(cy-sy)**two)**half          
          if((point<=(xn*yn)).and.(timea<=time<=timec))then
              if(timea<=time<=timeb)then
                  qp=pressure*(time/timeb)
              end if
              if(timeb<time<=timec)then
                  qp=pressure*(timec-time)/(timec-timeb)
              end if
          end if
          if((time>timec).or.(distance>r))then
              qp=0.d0
          end if
          value(km)=qp*exp(-(distance**two)/(two*r*r))    
      end do
      return
      end

​