/config,nproc,2
/FILNAME, dc_adapbounday
/PREP7
/title, directively coupled analysis with heat flux on top surface
! /OUTPUT,,output

ET,1,SOLID5
KEYOPT,1,1,0
MP,DENS,1,2700
MPTEMP,,0,50,100,150,200,250,300
MPDATA,KXX,1,,90,100,110,120,130,140,150
MPDATA,C,1,,850,900,950,1000,1050,1100,1120
mp,ex,1, 69e9
mp,alpx,1,0.0000236
mp,nuxy,1,.33 
TB,BKIN,1,,2,1
TBDATA,,241e6,0,,,,
whinput=1860
length=0.304
width=0.102
thick=0.0081
rshoulder=0.0127
hbackplate=0.002
nthke=4
smesize=0.004
ncopy=length/(4*smesize)
block,,4*smesize,,thick,,5*smesize 
block,,4*smesize,,thick,5*smesize,7*smesize
block,,4*smesize,,thick,7*smesize,9*smesize
block,,4*smesize,,thick,9*smesize,width
nummrg,all
numcmp,all
lesize,1,,,nthke
MSHAPE,0,3D
MSHKEY,1
esize,smesize 
vMESH,1

lesize,15,,,2,
lesize,13,,,2,
lesize,14,,,4,
lesize,16,,,4,
lesize,17,,,1,
lesize,20,,,1,
lesize,18,,,1,
lesize,19,,,1,
vMESH,2
lesize,22,,,2,
lesize,24,,,2,
lesize,21,,,2,
lesize,23,,,2, 
lesize,25,,,1,
lesize,26,,,1,
lesize,27,,,1,
lesize,28,,,1,
vmesh,3
lesize,36,,,8
vmesh,4
VGEN,ncopy, all,,,4*smesize,,,,,
nummrg,all
VSYMM, Z, all
nummrg,all
WPOFFS, (rshoulder+length/8), thick 
WPROTA,,-90
CYLIND,rshoulder, ,,2*thick,0,360,
WPOFFS, -(rshoulder+length/8), -thick
WPROTA,,90
MP,MU,1,0.3
MAT,1
R,3
REAL,3
ET,2,170
ET,3,173
KEYOPT,3,9,0
KEYOPT,3,10,2
R,3,
RMORE, 
RMORE,,0
RMORE,0
! Generate the target surface
asel,s,,,26
CM,_TARGET,AREA
AATT,-1,3,2,-1
TYPE,2
AMESH,ALL
allsel
! Create a pilot node
N,10000, (rshoulder+length/8),thick,0
TSHAP,PILO 
E,10000
! Generate the contact surface
asel,s,loc,y,thick
asel,r,loc,z,-5*smesize,5*smesize
asel,u,,,26
CM,_CONTACT,AREA
TYPE,3
NSLA,S,1
ESLN,S,0
ESURF
ALLSEL
block,-4*smesize,length+4*smesize,-thick-smesize,-thick,-width-4*smesize,width+4*smesize 
/COM, CONTACT PAIR CREATION - START
MP,MU,1,0.3
MAT,1
R,4
REAL,4
ET,4,170
ET,5,173
KEYOPT,5,1,1
KEYOPT,5,9,0
KEYOPT,5,10,2
R,4,
! Generate the target surface
ASEL,S,,,451
CM,_TARGET,AREA
AATT,-1,4,4,-1 
TYPE,4
AMESH,ALL
! Create a pilot node
N,16000, 0,0,0
TSHAP,PILO
E,16000
! Generate the contact surface
allsel
ASEL,S,loc,y,0
ASEL,U,,,451
CM,_CONTACT,AREA
TYPE,5
NSLA,S,1
ESLN,S,0 
ESURF
*SET,_REALID,4
ALLSEL
ESEL,ALL
/PNUM,TYPE,1
/NUM,1
EPLOT
/COM, CONTACT PAIR CREATION - END
FINISH
/SOLU
tref,25
allsel
esel,s,type,,1
nsle
nsel,r,loc,z,width
nsel,r,loc,y,0 
d,all,all,0
allsel
esel,s,type,,1
nsle
nsel,r,loc,z,-width
nsel,r,loc,y,0
d,all,all,0
allsel
esel,s,type,,1
nsle
nsel,r,loc,y,0
SF,all,CONV,350,20
allsel
esel,s,type,,1
nsle
nsel,r,loc,y,thick
SF,all,CONV,30,20
allsel
esel,s,type,,1
nsle
nsel,r,loc,x,0
SF,all,CONV,30,20
allsel
esel,s,type,,1
nsle
nsel,r,loc,x,length
SF,all,CONV,30,20
allsel
esel,s,type,,1 
nsle
nsel,r,loc,z,width
SF,all,CONV,30,20
allsel
esel,s,type,,1
nsle
nsel,r,loc,z,-width
SF,all,CONV,30,20
allsel
D,10000,UZ,0
D,10000,ROTX,0
D,10000,ROTZ,0
F,10000,FY,-10000
d,16000,all,0 
allsel
Antype,trans
NLGEOM,off
NROPT,UNSYM
KBC,0
OUTRES,,last
dtime=0.1                      !time step
vel=0.00336                       !welding speed
xinc=dtime*vel                  !increment in X direction
nloct=NINT(length/(80*xinc)) 
   esel,s,type,,1
   nsle
   nsel,r,loc,y,thick
   nsel,r,loc,z,-(rshoulder+smesize),(rshoulder+smesize)
   *get,numnode,node,0,count   !number of nodes in the selected set
   *get,lnnum,node,0,num,min    !loest node number in the selected set 
   *dim,nodenumber,ARRAY,numnode
   *dim,nodehflux,ARRAY,numnode
   *dim,emnumber,ARRAY,numnode
   allsel
      nsel,s,loc,y,0
      esln,s
      esel,r,type,,5
  *get,numelemback,ELEM,0,count   !number of elements in the selected set 
   *get,lenumback,ELEM,0,num,min    !loest element number in the selected set
   *dim,elemnumberback,ARRAY,numelemback
   *dim,contgap,ARRAY,numelemback
   *dim,elemconv,ARRAY,numelemback
   *dim,contpressure,ARRAY,numelemback
    backelemnumber=lenumback
   *do,ja,1,numelemback
      elemnumberback(ja)= backelemnumber
      backnoden=ELNEXT(backelemnumber)
      backelemnumber=backnoden
   *enddo
allsel
time,0.01
solve
save
finish
/solu
*do,i,1,nloct,1
   x0=rshoulder+(i-1)*xinc+length/8   !start point
   z0=0
   esel,s,type,,1
   nsle
   nsel,r,loc,y,thick
   nsel,r,loc,z,-(rshoulder+smesize),(rshoulder+smesize)
    nnum=lnnum
  	  jj=0
   	*do,ia,1,numnode
	 	*get,xlocation,node,nnum,loc,x
    	 	 *get,zlocation,node,nnum,loc,z
     		  rd=SQRT((xlocation-x0)**2+(zlocation-z0)**2)
       	 		 *if,rd,le,rshoulder,THEN
        	 			  jj=jj+1 
         				  nodenumber(jj)=nnum
        	 			  nodehflux(jj)=3*whinput*rd/(2*3.1416*rshoulder**3)
        			*endif
 	    	 *get,nexthn,node,nnum,nxth
   	 	   nnum=nexthn
  	  *enddo 
   allsel
             *do,ib,1,jj,1 
                     *if,ib,eq,1,then
                              nsel,s,node,,nodenumber(ib)
                     *else
                              nsel,a,node,,nodenumber(ib)
                              n3=nodenumber(ib)
                      *endif
              *enddo  
    ESLN,s,
    esel,r,type,,1
    *get,numelem,elem,0,count
    *get,lenum,elem,0,num,min
    enum=lenum
    nsle,all
           *do,ii,1,numelem
                 emnumber(ii)=enum
                 nnumk=nelem(enum,3)
                 nnuml=nelem(enum,4)
                 nnumo=nelem(enum,7) 
                 nnump=nelem(enum,8)
                 hfluxl=0
                 hfluxk=0
                 hfluxp=0
                 hfluxo=0
                       *do,iii,1,jj
                               *if,nnuml,eq,nodenumber(iii),then
                                       hfluxl=nodehflux(iii)
                               *elseif,nnumk,eq,nodenumber(iii) 
                                       hfluxk=nodehflux(iii)
                               *elseif,nnump,eq,nodenumber(iii)
                                       hfluxp=nodehflux(iii)
                              *elseif,nnumo,eq,nodenumber(iii)
                                       hfluxo=nodehflux(iii)
                              *endif 
                       *enddo 
                 aehflux=(hfluxl+hfluxk+hfluxp+hfluxo)/4
                 SFEDELE,enum,4,CONV
                 sfe,enum,4,HFLUX,,hfluxk,hfluxl,hfluxp,hfluxo
                 *get,nexthe,elem,enum,nxth
                 enum=nexthe
          *enddo 
 allsel 
  antype,trans,rest,i
   D,10000,UX,i*xinc,
   D,10000,ROTY,50*dtime*i,
   time,i*dtime+0.01
   DELTIM,dtime/5,dtime/50,dtime
   solve
   save
       /post1
      nsel,s,loc,y,0
      esln,s
      esel,r,type,,5
      ETABLE,ctpres,CONT,PRES
      ETABLE,ctgap,CONT,GAP
              *do,jb,1,numelemback
                   *GET, contpressure(jb), ctpres, elemnumberback(jb)
                   *GET, contgap(jb), ctgap, elemnumberback(jb)
                  *if,contpressure(jb)*0.000145,LE,100,and,contpressure(jb)*0.000145,GT,0,then
                        elemconv(jb)=(3600-3000)*(contpressure(jb)*0.000145-5)/(100-5)+3000
                  *elseif,contpressure(jb),LT,100
                        elemconv(jb)=(5000-3600)*(contpressure(jb)*0.000145-100)/(380-100)+3600
                  *else
                        elemconv(jb)=100      
                  *endif
              *enddo 
allsel
/solu
             *do,id,1,numelem
                     enum=emnumber(id)
                     SFEDELE,enum,4,hflux
             *enddo  
    allsel
    esel,s,type,,1
    nsle
    nsel,r,loc,y,thick
   SF,all,CONV,30,20
   allsel
            *do,jc,1,numelemback
                   sfe,elemnumberback(jc),1,CONV,,elemconv(jc)
                   sfe,elemnumberback(jc),1,CONV,2,25
            *enddo 
 allsel 
*enddo

F,10000,FY,0
D,10000,UX,i*xinc,
D,10000,UY,10*xinc,
nsel,all
nsel,s,loc,z,width
nsel,a,loc,z,-width
nsel,r,loc,y,0
DDELE,all,all
allsel
  antype,trans,rest,i+1
       time,i*dtime+0.02
   DELTIM,dtime/10,dtime/50,dtime
       solve
*do,je,1,10
   finish
    /post1
        nsel,s,loc,y,0
      esln,s
      esel,r,type,,5
      ETABLE,ctpres,CONT,PRES
      ETABLE,ctgap,CONT,GAP
   *do,jf,1,numelemback
    *GET, contpressure(jf), ctpres, elemnumberback(jf)
    *GET, contgap(jf), ctgap, elemnumberback(jf)
   *if,contpressure(jf)*0.000145,LE,100,and,contpressure(jf)*0.000145,GT,0,then
             elemconv(jf)=(3600-3000)*(contpressure(jf)*0.000145-5)/(100-5)+3000
     *elseif,contpressure(jf),LT,100
             elemconv(jf)=(5000-3600)*(contpressure(jf)*0.000145-100)/(380-100)+3600
     *else
            elemconv(jf)=100      
   *endif
  *enddo 
finish
/solu
   allsel
   *do,jg,1,numelemback
      sfe,elemnumberback(jg),1,CONV,,elemconv(jg)
      sfe,elemnumberback(jg),1,CONV,2,25
  *enddo 
 allsel 
  antype,trans,rest,i+1+je
       time,i*dtime+0.5*je+0.02
   DELTIM,dtime/10,dtime/50,dtime
       solve 
*enddo
finish