/BATCH  
/CONFIG,nproc,3
rn=103
/filnam,global2_aa_%rn%
!resume
/prep7
!*if,1,eq,1,:2000
!/input,test,cdb      ! Original model definition
!/input,kf3coil010407,cdb
!/input,kfibbolts020907,cdb    ! all possible inboard bolts
!/input,kfibbolts021507,cdb    ! PPPL choice inboard bolts
!/input,kfibbolts031307,cdb    ! PPPL choice inboard bolts
/input,kfibbolts032607,cdb    ! PPPL choice inboard bolts

allsel 
lsclear,all

/com
/com Global model of NCSX MCs from ORNL
/com Model Bolted flanges one at a time
/com Set k_bolts to one of the following:
/com   
/com
k_bolts=1    ! 1: AA, 2: AB, 3: BC, 4: CC
k_cont=4     ! 0: std, 1: rough, 4: sliding, 5: bonded
mu_cont=0.3  ! contact surface friction coef.

inner_friction = 1.3	!set inner leg friction

/com Number of bolts at each interface (comes from nodal components in cdb file)
n1_bolts=20  ! AA bolts
n2_bolts=29  ! A & B bolts
n3_bolts=33  ! BB & C bolts
n4_bolts=32  ! C bolts

/com
/com Run History
/com
/com 0: Add bolts, but keep flanges stuck (for run-time benchmark)
/com 1: Add A-B bolts, mu=0.3, preload from reft in bolts. Poor convergence
/com 2: A-B bolts, bonded, t(bolt)= -80, tref/tunif=300/300, PCG. Looks OK, but preload too light at 56 kips/bolt.
/com 3: A-B bolts, mu=0.00, t(bolt)=-200, tref/tunif=300/300, PCG, 50 substeps with autots. Preload converged nicely. EM died at 1.04 s.
/com 4: A-B bolts, mu=0.25, kfn=1, t(bolt)=-200, tref/tunif=300/300, PCG, 50 & 100 substeps with autots. 4.5 days 117 CUM, Preload + 4%of EM
/com 5: A-B bolts, mu=0.00, keyopt12=4 (sliding), kfn=1, t(bolt)=-200, tref/tunif=300/300, PCG, 50 & 100 substeps with autots, dbolt=1.15".
/com 6: A-B bolts, mu=0.00, keyopt12=4 (sliding), kfn=1, t(bolt)=-200, tref/tunif=300/300, PCG, 50 & 100 substeps with autots, dbolt=2.75".
/com 7: A-B bolts, mu=0.15, keyopt12=4 (sliding), kfn=1, t(bolt)=-200, tref/tunif=300/300, PCG, 50 & 100 substeps with autots, dbolt=2.75". 
/com     Would have converged in ~30 hr.

/com
/com Kevin Fixed Bolt Holes and made nodal components for all of them
/com  Make keyopt(12)=4 for sliding joint (with or without friction)
/com  Make preload with 500K deltaT
/com 210: mu=0.15 at A-B & BB-C: no conv. at Preload.
/com 210: mu=0.20 at A-B with corrected F(TF) loads. 2.75" with 500K Preload: ~1.25 MN (281k-lb). Too high. Should be 73000 lb.
/com 211: mu=0.20 at A-B with corrected F(TF) loads. 1.375" with 250K Preload  78 k-lb (17 hrs)
/com 221: mu=0.25 at A-B with corrected F(TF) loads. 1.375" with 250K Preload  78 k-lb (17 hrs)
/com 222: mu=0.20 in  bolt region of A-B, SS puck joint stiffness at inboard leg. 
/com BB-C
/com 212: mu=0.20 at BB-C with corrected F(TF) loads. 1.375" with 270K Preload 78 k-lb ave, light on end bolts (37 hrs)
/com 
/com A-A
/com 213: mu=0.0 at A-A & C-C, 1.375" with 300K Preload ~78 k-lb ave (2 hrs)
/com 213: mu=0.2 at A-A      , 1.375" with 300K Preload ~78 k-lb ave ( hrs)
/com 
/com Bond all interfaces for estimate of inboard leg shear
/com
/com 20: Bonded
/com

/com
/com Add possible inboard Bolts
/com
/com 230: AB, 0.25 friction, keyopt12=4
/com 231: AB, 0.00 friction, keyopt2/12=0/1 (hi-mu, but general contact): Poor convergence
/com 231: AB, 0.00 friction, keyopt2/12=1/1 (hi-mu, but general contact): Poor convergence
/com
/com PPPL Choice of Inboard Bolts
/com
/com 241: AA, 0.3 friction, keyopt12=4 (sliding)
/com 242: AB, 0.3 friction, keyopt12=4
/com 243: BC, 0.3 friction, keyopt12=4
/com 244: CC, 0.3 friction, keyopt12=4
/com 245: AA, 0.3 friction, keyopt12=1 (rough), sparse, unsym, update every iter: 3-4 iterations in 15 hrs. Killed
/com 245: AA, 0.3 friction, keyopt12=1 (rough),    PCG,   sym, update elem every iter:  t=0.4, 74 iterations in 48 hrs. Killed
/com 245: AA, 0.3 friction, keyopt12=0 (std),PCG: 1e-7, unsym, update pair every iter, fkn=0.1: 15 hr produced 
/com 246: AB, 0.3 friction, keyopt12=1
/com 247: BC, 0.3 friction, keyopt12=1
/com 248: CC, 0.3 friction, keyopt12=1

/com Study Transverse stiffness of contact elements

/com Revised AA
/com 251: AA, 0.4 friction, keyopt12=4 (sliding), 6 Inboard Bolts, specify Kt=1e11: Slip<0.20mm, T(bolt)<10kip
/com 252: AA, 0.4 friction, keyopt12=4 (sliding), 6 Inboard Bolts, specify Kt=5e11: Slip<0.12mm, T(bolt)<5kip
/com 253: AA, 0.4 friction, keyopt12=4 (sliding), 6 Inboard Bolts, specify Kt=10e11: Slip<0.11mm, T(bolt)<4.5kip

/com Looks like 5e11 is a reasonable value

/com 254: AB, 0.4 fric, keyopt12=4: Slip<?mm, T(bolt)<?kip
/com 255: CC, 0.4 fric, keyopt12=4: Pre~65kip, Slip<0.17mm, T(bolt)<5kip (adjust dT to fix Preload but not rerun yet)
/com 256: BC

/com
/com Misc Parameters
/com
itm=0.0254   ! in. to m
*afun,rad
pi=acos(-1)
*afun,deg
/pnum,mat,1
/num,1
/edge,1,1

/com
/com Add a vertical constraint at the bottom ~theta=30 
/com
esel,s,type,,10
nsle
csys,1
nsel,r,loc,y,28,32
nsel,r,loc,z,-.05,.05
*get,zmn_grnd,node,,mnloc,x
nsel,r,loc,x,zmn_grnd
*get,n_grnd,node,,num,min
d,n_grnd,uz

/com
/com Material Properties
/com
alpha_struc=� 13e-6 � � � � � � � � �  ! Thermal expansion coefficient of Structure (should be 13e-6 for SS at 85K)


/com Shell Structures A,B,C as well as shims
*do,i,10,16
mp,murx,i,1
mp,�  ex,i,(22/0.145)*1e9
mp,alpx,i,0
*enddo

! Wing shims properties
*do,i,7,9   
mp,murx,i,1
mp,�  ex,i,(1/0.145)*1e9  !weak = .0001, strong = 1
mp,alpx,i,0
*enddo

/com Bolts used for washers of the Bolted Joints
mp,murx,100,1
mp,prxy,100,.3
mp,ex,100,1000*190e9
mp,alpx,100,0

/com Bolts used for A-B bolts
mp,murx,101,1
mp,prxy,101,.3
mp,ex,101,190e9
mp,alpx,101,alpha_struc

/com Bolts used for B-C bolts
mp,murx,102,1
mp,prxy,102,.3
mp,ex,102,190e9
mp,alpx,102,alpha_struc

/com Bolts used for AA bolts
mp,murx,103,1
mp,prxy,103,.3
mp,ex,103,190e9
mp,alpx,103,alpha_struc

/com Bolts used for CC bolts
mp,murx,104,1
mp,prxy,104,.3
mp,ex,104,190e9
mp,alpx,104,alpha_struc
� 
! mod coil properties
*do,j,1,6
mp,ex,j,(8.5/100/0.145)*1e9
mp,alpx,j,(400e-6/40)
*enddo

/com
/com Reflect the shims about Z=0 and numm, shims should be last two element types
/com
allsel
*get,nmx,node,,num,max
esel,s,type,,15
nsle
*get,nmnshms15,node,,num,min
*get,nmxshms15,node,,num,max
csys
nsym,z,nmx-nmnshms15+1,all
esym, ,nmx-nmnshms15+1,all
numm,node
allsel
*get,nmx,node,,num,max
esel,s,type,,16
nsle
*get,nmnshms16,node,,num,min
*get,nmxshms16,node,,num,max
csys
nsym,z,nmx-nmnshms16+1,all
esym, ,nmx-nmnshms16+1,all
numm,node

/com
/com Add Bonded Target & Contact Elements between Shims and Shell Flanges
/com
esel,s,type,,16
nsle
csys
nsel,r,loc,y,0.006349,0.006351
cm,nshim_th0,node
esel,s,type,,11
nsle
nsel,r,loc,y,0.006349,0.006351
cm,nshell_th0,node
local,11,,,,,60
esel,s,type,,15
nsle
nsel,r,loc,y,-0.006351,-0.006349
cm,nshim_th60,node
esel,s,type,,12
nsle
nsel,r,loc,y,-0.006351,-0.006349
cm,nshell_th60,node


/com AA
/COM, CONTACT PAIR CREATION - START 
CM,_NODECM,NODE 
CM,_ELEMCM,ELEM 
CM,_KPCM,KP 
CM,_LINECM,LINE 
CM,_AREACM,AREA 
CM,_VOLUCM,VOLU 
/GSAV,cwz,gsav,,temp
MP,MU,66,
MAT,66   
R,65
REAL,65 
ET,65,170   
ET,66,174   
R,65,,,1.0,0,0,   
R,65,,,0.1,0,0,   
!R,65,,,1.0,0.1,0,   
RMORE,,,1.0E20,0.0,1.0,-5e11          ! define shear stiffness since default too soft
!RMORE,0.0,0,1.0,,1.0,0.5
!RMORE,0,1.0,1.0,0.0,,1.0
KEYOPT,66,4,0   
KEYOPT,66,5,0   
KEYOPT,66,7,0   
KEYOPT,66,8,0   
KEYOPT,66,9,1   ! Exclude both initial geometrical penetration or gap and offset
KEYOPT,66,10,5 !1  
KEYOPT,66,11,0  
KEYOPT,66,12,5  
KEYOPT,66,2,0   
KEYOPT,65,5,0   
! Generate the target surface   
NSEL,S,,,NSHELL_TH0 
CM,_TARGET,NODE 
TYPE,65 
ESLN,S,0
ESURF   
CMSEL,S,_ELEMCM 
! Generate the contact surface  
NSEL,S,,,NSHIM_TH0  
CM,_CONTACT,NODE
TYPE,66 
ESLN,S,0
ESURF   
ALLSEL  
ESEL,ALL
ESEL,S,TYPE,,65 
ESEL,A,TYPE,,66 
ESEL,R,REAL,,65 
/PSYMB,ESYS,1   
/PNUM,TYPE,1
/NUM,1  
!EPLOT   
ESEL,ALL
ESEL,S,TYPE,,65 
ESEL,A,TYPE,,66 
ESEL,R,REAL,,65 
CMSEL,A,_NODECM 
CMDEL,_NODECM   
CMSEL,A,_ELEMCM 
CMDEL,_ELEMCM   
CMSEL,S,_KPCM   
CMDEL,_KPCM 
CMSEL,S,_LINECM 
CMDEL,_LINECM   
CMSEL,S,_AREACM 
CMDEL,_AREACM   
CMSEL,S,_VOLUCM 
CMDEL,_VOLUCM   
/GRES,cwz,gsav  
CMDEL,_TARGET   
CMDEL,_CONTACT  
/COM, CONTACT PAIR CREATION - END
 
/com CC 
/COM, CONTACT PAIR CREATION - START 
CM,_NODECM,NODE 
CM,_ELEMCM,ELEM 
CM,_KPCM,KP 
CM,_LINECM,LINE 
CM,_AREACM,AREA 
CM,_VOLUCM,VOLU 
/GSAV,cwz,gsav,,temp
MP,MU,64,0   
MAT,64
R,66
REAL,66 
ET,63,170   
ET,64,174   
R,66,,,1.0,0,0,   
R,66,,,0.1,0,0,   
RMORE,,,1.0E20,0.0,1.0,-5e11          ! define shear stiffness since default too soft
!RMORE,0.0,0,1.0,,1.0,0.5
!RMORE,0,1.0,1.0,0.0,,1.0
KEYOPT,64,4,0   
KEYOPT,64,5,0   
KEYOPT,64,7,0   
KEYOPT,64,8,0   
KEYOPT,64,9,1     ! Exclude both initial geometrical penetration or gap and offset
KEYOPT,64,10,5 !1  
KEYOPT,64,11,0  
KEYOPT,64,12,5  
KEYOPT,64,2,0   
KEYOPT,63,5,0   
! Generate the target surface   
NSEL,S,,,NSHELL_TH60
CM,_TARGET,NODE 
TYPE,63 
ESLN,S,0
ESURF   
CMSEL,S,_ELEMCM 
! Generate the contact surface  
NSEL,S,,,NSHIM_TH60 
CM,_CONTACT,NODE
TYPE,64 
ESLN,S,0
ESURF   
ALLSEL  
ESEL,ALL
ESEL,S,TYPE,,63 
ESEL,A,TYPE,,64 
ESEL,R,REAL,,66 
/PSYMB,ESYS,1   
/PNUM,TYPE,1
/NUM,1  
!EPLOT   
ESEL,ALL
ESEL,S,TYPE,,63 
ESEL,A,TYPE,,64 
ESEL,R,REAL,,66 
CMSEL,A,_NODECM 
CMDEL,_NODECM   
CMSEL,A,_ELEMCM 
CMDEL,_ELEMCM   
CMSEL,S,_KPCM   
CMDEL,_KPCM 
CMSEL,S,_LINECM 
CMDEL,_LINECM   
CMSEL,S,_AREACM 
CMDEL,_AREACM   
CMSEL,S,_VOLUCM 
CMDEL,_VOLUCM   
/GRES,cwz,gsav  
CMDEL,_TARGET   
CMDEL,_CONTACT  
/COM, CONTACT PAIR CREATION - END

/com remove the friggin shim part not supported

/com gets the inner leg no go region on coil A.

esel,s,type,,16
nsle

angle = 0

csys,0
wpcs
wpoffs,1.5
wprota,angle,90
/view,,wp

cswpla,1000,1
nsel,r,loc,y,135,158
esln
cm,upper_A,elem

esel,s,type,,16
nsle
nsel,r,loc,y,135,360-135
esln
cm,lower_A,elem
cmgrp,nobolts_A,lower_A,upper_A

esel,s,,,nobolts_A
nsle
esln
esel,r,type,,66
cm,slideAA1,elem


csys,0
wpcs
wpoffs,.8
wprota,angle,90
/view,,wp
esel,s,type,,16
nsle
cswpla,1000
nsel,r,loc,x,0,1
nsel,r,loc,y,-1,1
esln
esel,r,type,,66
cm,slideAA2,elem
cmgrp,slideaa,slideAA2,slideAA1


/com done with this, do for ab bc and cc if needed.


/com 
/com A-to-A CEs for Three-Coil Structural Model
/com
tol=0.0001� � � � �  ! tolerance for pairing nodes for up-down CEs
esel,s,mat,,16
nsle
csys
nsel,r,loc,y,-tol,tol
cm,n_asym,node
nsel,r,loc,z,0.0001,12
cm,npz_asym,node
*get,n_strt,node,,num,min
*get,n_stop,node,,num,max
*do,j,n_strt,n_stop,1 !5
cmsel,s,npz_asym
*if,nsel(j),ne,1,cycle
cmsel,s,n_asym
nsel,r,loc,x,nx(j)-tol,nx(j)+tol
nsel,r,loc,z,-nz(j)-tol,-nz(j)+tol
*get,any,node,,count
*if,any,eq,0,cycle
*get,j_asym,node,,num,min
ce,next,,j,ux,1,j_asym,ux,-1� �  ! Asymmetric nodes move in X together
ce,next,,j,uy,1,j_asym,uy,+1� �  ! Asymmetric nodes move in opposite Y
ce,next,,j,uz,1,j_asym,uz,+1� �  ! Asymmetric nodes move in opposite Z
*enddo

/com 
/com C-to-C CEs for Three-Coil Structural Model
/com
esel,s,mat,,15
nsle
local,11,,,,,60
nsel,r,loc,y,-tol,tol
cm,n_asym,node
nrotate,all
nsel,r,loc,z,.0001,12
cm,npz_asym,node
*get,n_strt,node,,num,min
*get,n_stop,node,,num,max
cmsel,s,n_asym
*do,j,n_strt,n_stop,1
cmsel,s,npz_asym
*if,nsel(j),ne,1,cycle
cmsel,s,n_asym
nsel,r,loc,x,nx(j)-tol,nx(j)+tol
nsel,r,loc,z,-nz(j)-tol,-nz(j)+tol
*get,any,node,,count
*if,any,eq,0,cycle
*get,j_asym,node,,num,min
ce,next,,j,ux,1,j_asym,ux,-1� �  ! Asymmetric nodes move in X together
ce,next,,j,uy,1,j_asym,uy,+1� �  ! Asymmetric nodes move in opposite Y
ce,next,,j,uz,1,j_asym,uz,+1� �  ! Asymmetric nodes move in opposite Z
*enddo
csys
allsel


/com APPLY TEMPERATURES (not neeeded for now.)
TREF,300
TUNIF,300
!/com esel,s,type,,2,7
!/com bfe,all,temp,,120
allsel


!!!!select out wings if applicable
!esel,u,type,,7,9	! wing ele's	
!esel,u,type,,35,48  	! contact ele's for wings

/com
/com Flanges pulled together by rods through holes (type 101, 102, 103 & 104 for AB, BC, AA & CC).
/com End of rods connected to hole perimeter by little pipes (type 100)
d_bolt=2.75*itm     ! when slip is high, large bolt has equivalent stiffness to collared bolt
d_bolt=1.375*itm    ! for predominantly tensile loaded bolts with friction and minimal slip
a_bolt=(pi/4)*d_bolt**2
i_bolt=(pi/64)*d_bolt**4
j_bolt=(pi/32)*d_bolt**4
/com Spokes
et,100,pipe16
r,100,d_bolt,d_bolt/2
/com A-B Joints
et,101,pipe16
r,101,d_bolt,d_bolt/2
/com BB-C Joints
et,102,pipe16
r,102,d_bolt,d_bolt/2
/com A-A Joints
et,103,pipe16
r,103,d_bolt,d_bolt/2
/com C-C Joints
et,104,pipe16
r,104,d_bolt,d_bolt/2

/com Nodes to indicate which way the local CSYS is pointing
nminusy=999998
nplusy=999999
csys
n,nminusy,,-2
n,nplusy,,2

/com Flange A/B 
A = 1000000
*do,j,1,n2_bolts
cmsel,s,holea%j%,node
*get,ncount,node,,count
*do,i,1,ncount
*get,n%i%,node,,num,min
nsel,u,,,n%i%
*enddo

cmsel,s,holea%j%,node
nwplan,,n1,n2,n3
*get,csmx,cdsy,,num,max
cswpla,csmx+1

!coords of points
n1x = NX(n1)
n1y = NY(n1)
n1z = NZ(n1)

*do,k,2,ncount
*if,NX(n1),ne,NX(n%k%),and,Ny(n1),ne,Ny(n%k%),exit
*enddo

n2x = NX(n%k%)
n2y = NY(n%k%)
n2z = Nz(n%k%)

*do,m,k,ncount
*if,n2x,ne,NX(n%m%),and,n2y,ne,Ny(n%m%),exit
*enddo

n3x = NX(n%m%)
n3y = NY(n%m%)
n3z = Nz(n%m%)

ma = (n2y-n1y)/(n2x-n1x)
mb = (n3y-n2y)/(n3x-n2x)

xc = (ma*mb*(n1y-n3y)+mb*(n1x+n2x)-ma*(n2x+n3x))/2/(mb-ma)
yc = -1/ma*(xc-(n1x+n2x)/2)+(n1y+n2y)/2
ycc = -1/mb*(xc-(n1x+n2x)/2)+(n1y+n2y)/2

/com make node at center
*if,ma,lt,1e-6,and,ma,gt,-1e-6,then
n,A,xc,ycc,n3z
*else
n,A,xc,yc,n3z
*endif
!cp,next,all,A,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12     ! Nodes at flange face near shim
/com Make WP at A
nwpla,,A,n2,n3
/com Check Orientation, Offset to Washer face & make csys
cswpla,csmx+1,1
*if,nz(nminusy),gt,0,then
wprot,,180
*endif
/com Offset to Washer face & make csys
wpoff,,,-1.25*0.0254
cswpla,csmx+1,1
/com Check to see if this is a viable fastener seat
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,a
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*if,nn,lt,12,then
wpoff,,,1.25*0.0254
cswpla,csmx+1,1
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,a
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*endif
cm,n_a%j%,node
n,A
/com Connect center node to perimeter nodes with little pipes
mat,100
type,100
real,100
*do,i,1,nn
*if,nn,lt,12,cycle      ! check for holes without well-defined nut-side surfaces
*get,nmin,node,,num,min
e,A,nmin
nsel,u,node,,nmin
*enddo
A = A + 1
*enddo
! end flange A

/com Flange B
B = 2000000
*do,j,1,n2_bolts
cmsel,s,holeb%j%,node
*get,ncount,node,,count
*do,i,1,ncount
*get,n%i%,node,,num,min
nsel,u,,,n%i%
*enddo

cmsel,s,holeb%j%,node
nwplan,,n1,n2,n3
*get,csmx,cdsy,,num,max
cswpla,csmx+1

!coords of points
n1x = NX(n1)
n1y = NY(n1)
n1z = Nz(n1)

*do,k,2,ncount
*if,NX(n1),ne,NX(n%k%),and,Ny(n1),ne,Ny(n%k%),exit
*enddo

n2x = NX(n%k%)
n2y = NY(n%k%)
n2z = Nz(n%k%)

*do,m,k,ncount
*if,n2x,ne,NX(n%m%),and,n2y,ne,Ny(n%m%),exit
*enddo

n3x = NX(n%m%)
n3y = NY(n%m%)
n3z = Nz(n%m%)

ma = (n2y-n1y)/(n2x-n1x)
mb = (n3y-n2y)/(n3x-n2x)

xc = (ma*mb*(n1y-n3y)+mb*(n1x+n2x)-ma*(n2x+n3x))/2/(mb-ma)
yc = -1/ma*(xc-(n1x+n2x)/2)+(n1y+n2y)/2
ycc = -1/mb*(xc-(n1x+n2x)/2)+(n1y+n2y)/2

/com make node at the center
*if,ma,lt,1e-6,and,ma,gt,-1e-6,then
n,B,xc,ycc,n3z
*else
n,B,xc,yc,n3z
*endif
!cp,next,all,B,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12
/com Make WP at B
nwpla,,B,n2,n3
/com Check Orientation, Offset to Washer face & make csys
cswpla,csmx+1,1
*if,nz(nplusy),gt,0,then
wprot,,180
*endif
/com Offset to Washer face & make csys
wpoff,,,-1.25*0.0254
cswpla,csmx+1,1
/com Check to see if this is a viable fastener seat
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,b
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*if,nn,lt,12,then
wpoff,,,1.25*0.0254
cswpla,csmx+1,1
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,b
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*endif
cm,n_b%j%,node
n,B
/com Connect center node to perimeter nodes with little pipes
mat,100
type,100
real,100
*do,i,1,nn
*if,nn,lt,12,cycle      ! check for holes without well-defined nut-side surfaces
*get,nmin,node,,num,min
e,B,nmin
nsel,u,node,,nmin
*enddo
B = B + 1
*enddo
/com
/com Add the Bolts between A & B
/com
allsel
*do,j,1,n2_bolts
mat,101
type,101
real,101
e,1000000+(j-1),2000000+(j-1)
*enddo

/com
/com Flange BB & C
/com
BB = 3000000
*do,j,1,n3_bolts
cmsel,s,holebb%j%,node
*get,ncount,node,,count
*do,i,1,ncount
*get,n%i%,node,,num,min
nsel,u,,,n%i%
*enddo

cmsel,s,holebb%j%,node
nwplan,,n1,n2,n3
*get,csmx,cdsy,,num,max
cswpla,csmx+1

!coords of points
n1x = NX(n1)
n1y = NY(n1)
n1z = NZ(n1)

*do,k,2,ncount
*if,NX(n1),ne,NX(n%k%),and,Ny(n1),ne,Ny(n%k%),exit
*enddo

n2x = NX(n%k%)
n2y = NY(n%k%)
n2z = Nz(n%k%)

*do,m,k,ncount
*if,n2x,ne,NX(n%m%),and,n2y,ne,Ny(n%m%),exit
*enddo

n3x = NX(n%m%)
n3y = NY(n%m%)
n3z = Nz(n%m%)

ma = (n2y-n1y)/(n2x-n1x)
mb = (n3y-n2y)/(n3x-n2x)

xc = (ma*mb*(n1y-n3y)+mb*(n1x+n2x)-ma*(n2x+n3x))/2/(mb-ma)
yc = -1/ma*(xc-(n1x+n2x)/2)+(n1y+n2y)/2
ycc = -1/mb*(xc-(n1x+n2x)/2)+(n1y+n2y)/2

/com make node at center
*if,ma,lt,1e-6,and,ma,gt,-1e-6,then
n,BB,xc,ycc,n3z
*else
n,BB,xc,yc,n3z
*endif
!cp,next,all,BB,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12     ! Nodes at flange face near shim
/com Make WP at BB
nwpla,,BB,n2,n3
/com Check Orientation, Offset to Washer face & make csys
cswpla,csmx+1,1
*if,nz(nminusy),gt,0,then
wprot,,180
*endif
/com Offset to Washer face & make csys
wpoff,,,-1.25*0.0254
cswpla,csmx+1,1
/com Check to see if this is a viable fastener seat
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,BB
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*if,nn,lt,12,then
wpoff,,,1.25*0.0254
cswpla,csmx+1,1
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,BB
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*endif
cm,n_bb%j%,node
n,BB
/com Connect center node to perimeter nodes with little pipes
mat,100
type,100
real,100
*do,i,1,nn
*if,nn,lt,12,cycle      ! check for holes without well-defined nut-side surfaces
*get,nmin,node,,num,min
e,BB,nmin
nsel,u,node,,nmin
*enddo
BB = BB + 1
*enddo
! end flange BB

/com Flange C
C = 4000000
*do,j,1,n3_bolts
cmsel,s,holec%j%,node
*get,ncount,node,,count
*do,i,1,ncount
*get,n%i%,node,,num,min
nsel,u,,,n%i%
*enddo

cmsel,s,holec%j%,node
nwplan,,n1,n2,n3
*get,csmx,cdsy,,num,max
cswpla,csmx+1

!coords of points
n1x = NX(n1)
n1y = NY(n1)
n1z = Nz(n1)

*do,k,2,ncount
*if,NX(n1),ne,NX(n%k%),and,Ny(n1),ne,Ny(n%k%),exit
*enddo

n2x = NX(n%k%)
n2y = NY(n%k%)
n2z = Nz(n%k%)

*do,m,k,ncount
*if,n2x,ne,NX(n%m%),and,n2y,ne,Ny(n%m%),exit
*enddo

n3x = NX(n%m%)
n3y = NY(n%m%)
n3z = Nz(n%m%)

ma = (n2y-n1y)/(n2x-n1x)
mb = (n3y-n2y)/(n3x-n2x)

xc = (ma*mb*(n1y-n3y)+mb*(n1x+n2x)-ma*(n2x+n3x))/2/(mb-ma)
yc = -1/ma*(xc-(n1x+n2x)/2)+(n1y+n2y)/2
ycc = -1/mb*(xc-(n1x+n2x)/2)+(n1y+n2y)/2

/com make node at the center
*if,ma,lt,1e-6,and,ma,gt,-1e-6,then
n,C,xc,ycc,n3z
*else
n,C,xc,yc,n3z
*endif
!cp,next,all,C,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12
/com Make WP at C
nwpla,,C,n2,n3
/com Check Orientation, Offset to Washer face & make csys
cswpla,csmx+1,1
*if,nz(nplusy),gt,0,then
wprot,,180
*endif
/com Offset to Washer face & make csys
wpoff,,,-1.25*0.0254
cswpla,csmx+1,1
/com Check to see if this is a viable fastener seat
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,b
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*if,nn,lt,12,then
wpoff,,,1.25*0.0254
cswpla,csmx+1,1
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,b
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*endif
cm,n_c%j%,node
n,C

/com Connect center node to perimeter nodes with little pipes
mat,100
type,100
real,100
*do,i,1,nn
*if,nn,lt,12,cycle      ! check for holes without well-defined nut-side surfaces
*get,nmin,node,,num,min
e,C,nmin
nsel,u,node,,nmin
*enddo
C = C + 1
*enddo

/com
/com Add the Bolts between BB & C
/com
allsel
*do,j,1,n3_bolts
mat,102
type,102
real,102
e,3000000+(j-1),4000000+(j-1)
*enddo

/com
/com Flange AA
/com
csys
AA = 5000000
*do,j,1,n1_bolts
cmsel,s,holeaa%j%,node
*get,ncount,node,,count
*do,i,1,ncount
*get,n%i%,node,,num,min
nsel,u,,,n%i%
*enddo

cmsel,s,holeaa%j%,node
nwplan,,n1,n2,n3
*get,csmx,cdsy,,num,max
cswpla,csmx+1

!coords of points
n1x = NX(n1)
n1y = NY(n1)
n1z = NZ(n1)

*do,k,2,ncount
*if,NX(n1),ne,NX(n%k%),and,Ny(n1),ne,Ny(n%k%),exit
*enddo

n2x = NX(n%k%)
n2y = NY(n%k%)
n2z = Nz(n%k%)

*do,m,k,ncount
*if,n2x,ne,NX(n%m%),and,n2y,ne,Ny(n%m%),exit
*enddo

n3x = NX(n%m%)
n3y = NY(n%m%)
n3z = Nz(n%m%)

ma = (n2y-n1y)/(n2x-n1x)
mb = (n3y-n2y)/(n3x-n2x)

xc = (ma*mb*(n1y-n3y)+mb*(n1x+n2x)-ma*(n2x+n3x))/2/(mb-ma)
yc = -1/ma*(xc-(n1x+n2x)/2)+(n1y+n2y)/2
ycc = -1/mb*(xc-(n1x+n2x)/2)+(n1y+n2y)/2

/com make node at center
*if,ma,lt,1e-6,and,ma,gt,-1e-6,then
n,AA,xc,ycc,n3z
*else
n,AA,xc,yc,n3z
*endif
!cp,next,all,AA,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12     ! Nodes at flange face near shim
/com Make WP at AA
nwpla,,AA,n2,n3
/com Check Orientation, Offset to Washer face & make csys
cswpla,csmx+1,1
*if,nz(nminusy),gt,0,then
wprot,,180
*endif
wpoff,,,1.25*0.0254
cswpla,csmx+1,1
/com Check to see if this is a viable fastener seat
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,AA
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*if,nn,lt,12,then
wpoff,,,-1.25*0.0254
cswpla,csmx+1,1
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,AA
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*endif
cm,n_aa%j%,node
n,AA
/com Then make a corresponding node at symmetry plane
csys,1
n,AA+100,nx(AA),0,nz(AA)
/com Connect center node to perimeter nodes with little pipes
mat,100
type,100
real,100
*do,i,1,nn
*if,nn,lt,12,cycle      ! check for holes without well-defined nut-side surfaces
*get,nmin,node,,num,min
e,AA,nmin
nsel,u,node,,nmin
*enddo
mat,103
type,103
real,103
e,AA,AA+100
AA = AA + 1
*enddo


/com
/com Flange CC
/com
CC = 6000000
*do,j,1,n4_bolts
cmsel,s,holecc%j%,node
*get,ncount,node,,count
*do,i,1,ncount
*get,n%i%,node,,num,min
nsel,u,,,n%i%
*enddo

cmsel,s,holecc%j%,node
nwplan,,n1,n2,n3
*get,csmx,cdsy,,num,max
cswpla,csmx+1

!coords of points
n1x = NX(n1)
n1y = NY(n1)
n1z = NZ(n1)

*do,k,2,ncount
*if,NX(n1),ne,NX(n%k%),and,Ny(n1),ne,Ny(n%k%),exit
*enddo

n2x = NX(n%k%)
n2y = NY(n%k%)
n2z = Nz(n%k%)

*do,m,k,ncount
*if,n2x,ne,NX(n%m%),and,n2y,ne,Ny(n%m%),exit
*enddo

n3x = NX(n%m%)
n3y = NY(n%m%)
n3z = Nz(n%m%)

ma = (n2y-n1y)/(n2x-n1x)
mb = (n3y-n2y)/(n3x-n2x)

xc = (ma*mb*(n1y-n3y)+mb*(n1x+n2x)-ma*(n2x+n3x))/2/(mb-ma)
yc = -1/ma*(xc-(n1x+n2x)/2)+(n1y+n2y)/2
ycc = -1/mb*(xc-(n1x+n2x)/2)+(n1y+n2y)/2

/com make node at center
*if,ma,lt,1e-6,and,ma,gt,-1e-6,then
n,CC,xc,ycc,n3z
*else
n,CC,xc,yc,n3z
*endif
!cp,next,all,CC,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12     ! Nodes at flange face near shim
/com Make WP at CC
nwpla,,CC,n2,n3
/com Check Orientation, Offset to Washer face & make csys
cswpla,csmx+1,1
*if,nz(nplusy),gt,0,then
wprot,,180
*endif
wpoff,,,1.25*0.0254
cswpla,csmx+1,1
/com Check to see if this is a viable fastener seat
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,CC
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*if,nn,lt,12,then
wpoff,,,-1.25*0.0254
cswpla,csmx+1,1
esel,s,type,,10,12
nsle
nsel,r,loc,z,-0.001,0.001
nsel,u,node,,CC
*get,rmn,node,,mnloc,x
nsel,r,loc,x,,1.05*rmn
*get,nn,node,,count
*endif
cm,n_cc%j%,node
n,CC
/com Then make a corresponding node at symmetry plane
csys,1
n,CC+100,nx(CC),60,nz(CC)
/com Connect center node to perimeter nodes with little pipes
mat,100
type,100
real,100
*do,i,1,nn
*if,nn,lt,12,cycle      ! check for holes without well-defined nut-side surfaces
*get,nmin,node,,num,min
e,CC,nmin
nsel,u,node,,nmin
*enddo
mat,104
type,104
real,104
e,CC,CC+100
CC = CC + 1
*enddo

/com 
/com A-to-A CEs for Bolts
/com
tol=0.0015� � � � �  ! tolerance for pairing nodes for up-down CEs
esel,s,mat,,103
nsle
csys
nsel,r,loc,y,-tol,tol
cm,n_asym,node
nsel,r,loc,z,0.0001,12
cm,npz_asym,node
*get,n_strt,node,,num,min
*get,n_stop,node,,num,max
*do,j,n_strt,n_stop,1
cmsel,s,npz_asym
*if,nsel(j),ne,1,cycle
cmsel,s,n_asym
nsel,r,loc,x,nx(j)-tol,nx(j)+tol
nsel,r,loc,z,-nz(j)-tol,-nz(j)+tol
*get,any,node,,count
*if,any,eq,0,cycle
*get,j_asym,node,,num,min
ce,next,,j,ux,1,j_asym,ux,-1� �  ! Asymmetric nodes move in X together
ce,next,,j,uy,1,j_asym,uy,+1� �  ! Asymmetric nodes move in opposite Y
ce,next,,j,uz,1,j_asym,uz,+1� �  ! Asymmetric nodes move in opposite Z
*enddo
/com Symmetry plane bolt has UY=0
esel,s,type,,103
nsle
nsel,r,loc,y
nsel,r,loc,z
!d,all,uy
� 
/com 
/com C-to-C CEs for Bolts
/com
esel,s,type,,104
nsle
local,11,,,,,60
nsel,r,loc,y,-tol,tol
cm,n_asym,node
nrotate,all
nsel,r,loc,z,.0001,12
cm,npz_asym,node
*get,n_strt,node,,num,min
*get,n_stop,node,,num,max
cmsel,s,n_asym
*do,j,n_strt,n_stop
cmsel,s,npz_asym
*if,nsel(j),ne,1,cycle
cmsel,s,n_asym
nsel,r,loc,x,nx(j)-tol,nx(j)+tol
nsel,r,loc,z,-nz(j)-tol,-nz(j)+tol
*get,any,node,,count
*if,any,eq,0,cycle
*get,j_asym,node,,num,min
ce,next,,j,ux,1,j_asym,ux,-1� �  ! Asymmetric nodes move in X together
ce,next,,j,uy,1,j_asym,uy,+1� �  ! Asymmetric nodes move in opposite Y
ce,next,,j,uz,1,j_asym,uz,+1� �  ! Asymmetric nodes move in opposite Z
*enddo

/com Nix types 57 & 58 (bottom of AA)...seem to be redundant to 65/66
esel,s,type,,57,58
edele,all
/com Nix types 59 & 60 (bottom of CC)...seem to be redundant to 63/64
esel,s,type,,59,60
edele,all

/com
/com AB Inboard Contact: type149,real49
/com
!et,149,174,,1
!KEYOPT,149,9,1  
!KEYOPT,149,12,5 
!R,149,,,, ,,
!RMORE,,,, ,,-2.7e11 
!mp,mu,149,0
!phiabmn=145   ! angular start of inboard leg
!phiabmx=215   ! angular   end of inboard leg
!esel,s,real,,49
!nsle
!local,2003,3,,,,,,,2
!nsel,r,loc,z,phiabmn,phiabmx
!esln
!esel,r,real,,49
!cm,e_ibab,elem
!emodif,all,real,149
!emodif,all,mat,149
!esel,r,type,,49
!emodif,all,type,149

/com
/com Change Selected Flange Contact elements from Bonded to Sliding, set mu, and preload bolts
/com

*if,k_bolts,eq,1,then
/com A-A
esel,s,type,,103
BFE,all,TEMP,1,-27,-27     ! WRT 300K reference temp
nsel,s,node,,5000120,5000122
nsel,a,node,,5000110
esln
BFE,all,TEMP,1,-710,-710     ! WRT 300K reference temp
/com Since Bolts at all flanges are present, Nix Bolts not required for analysis
esel,s,type,,101,102
esel,a,type,,104
nsle
esln
esel,r,type,,100,104
nsle
nsel,r,node,,1000000,10000000
edele,all
ndele,all
/com Set Proper contact element keyoptions
keyopt,66,2,0          ! 0 default 1 Penalty Function
!KEYOPT,66,5,3          ! Close gap/reduce penetration with auto CNOF
!KEYOPT,66,7,1          ! Automatic bisection of increment
KEYOPT,66,9,1          ! Include both initial geometrical penetration or gap and offset
keyopt,66,10,5         ! Each iteration based on current mean stress of underlying elements (pair based)
keyopt,66,12,k_cont    ! 0 (normal) 1 (rough) 4 (sliding), 5 (bonded)
mp,mu,66,mu_cont
mp,mu,1000,inner_friction
esel,s,,,slideaa
mpchg,1000,all
*endif

*if,k_bolts,eq,2,then
/com A-B
esel,s,type,,101
BFE,all,TEMP,1,-100,-100   ! WRT 300K reference temp
nsel,s,node,,1000005
nsel,a,node,,1000010,1000011
nsel,a,node,,1000016
nsel,a,node,,1000026,1000027
esln
esel,r,type,,101
BFE,all,TEMP,1,-200,-200     ! WRT 300K reference temp

/com Since Bolts at all flanges are present, Nix Bolts not required for analysis
esel,s,type,,102,104
nsle
esln
esel,r,type,,100,104
nsle
nsel,r,node,,1000000,10000000
edele,all
ndele,all
/com Set Proper contact element keyoptions
keyopt,49,2,0          ! 0 default 1 Penalty Function
KEYOPT,49,9,1          ! Exclude both initial geometrical penetration or gap and offset
keyopt,49,10,5         ! Update contact stiffness for each element after each iteration
keyopt,49,12,k_cont    ! 0 (normal) 1 (rough) 4 (sliding), 5 (bonded)
mp,mu,49,mu_cont
/com Change A-B Shear stiffness
R,49,,,.1,,,
RMORE,,,,,,-5e11     ! effective shear stiffness of interface
R,50,,,.1,,,
RMORE,,,,,,-5e11     ! effective shear stiffness of interface
*endif

*if,k_bolts,eq,3,then
/com BB-C
esel,s,type,,102
BFE,all,TEMP,1,-100,-100   ! WRT 300K reference temp
/com Since Bolts at all flanges are present, Nix Bolts not required for analysis
esel,s,type,,101
esel,a,type,,103,104
nsle
esln
esel,r,type,,100,104
nsle
nsel,r,node,,1000000,10000000
edele,all
ndele,all
/com Set Proper contact element keyoptions
keyopt,55,2,0          ! 0 default 1 Penalty Function
KEYOPT,55,9,1          ! Exclude both initial geometrical penetration or gap and offset
keyopt,55,10,5         ! Update contact stiffness for each element after each iteration
keyopt,55,12,k_cont    ! 0 (normal) 4 (sliding), 5 (bonded)
mp,mu,55,mu_cont
/com Change B-C Shear stiffness
R,55,,,.1,,,
RMORE,,,,,,-5e11     ! effective shear stiffness of interface
R,56,,,.1,,,
RMORE,,,,,,-5e11     ! effective shear stiffness of interface

*endif

*if,k_bolts,eq,4,then
/com C-C
esel,s,type,,104
BFE,all,TEMP,1,-30,-30       ! WRT 300K reference temp
nsel,s,node,,6000000,6000002
nsel,a,node,,6000029,6000031
esln
esel,r,type,,104
BFE,all,TEMP,1,-210,-210     ! WRT 300K reference temp
nsel,s,node,,6000015,6000016
esln
esel,r,type,,104
BFE,all,TEMP,1,-190,-190   ! WRT 300K reference temp

/com Since Bolts at all flanges are present, Nix Bolts not required for analysis
esel,s,type,,101,103
nsle
esln
esel,r,type,,100,104
nsle
nsel,r,node,,1000000,10000000
edele,all
ndele,all
/com Set Proper contact element keyoptions
keyopt,64,2,0          ! 0 default 1 Penalty Function
KEYOPT,64,9,1          ! Exclude both initial geometrical penetration or gap and offset
keyopt,64,10,5         ! Update contact stiffness for each element after each iteration
keyopt,64,12,k_cont    ! 0 (normal) 4 (sliding), 5 (bonded)
mp,mu,64,mu_cont
*endif

allsel
save
fini

:1000
/SOLU  
antyp
EQSLV,pcg,1e-6
!EQSLV,sparse
autots,on
nsubst,5,10,2
cnvtol,f,1000,0.1
!nropt,unsym    ! with friction
!MSAVE, on
!PRECISION, 1
!/com BCSOPTION,,,,,,PERFORMANCE
allsel

!esel,u,num,,314079
!esel,u,num,,314096
!esel,u,num,,314113
!esel,u,num,,314130
!esel,u,num,,314147
!esel,u,num,,314164

SOLVE

!:2000
!/solu  
!ANTYPE,,REST, , ,0  
! pulls in mod coil em forces solved by seperate em model provided by Len
/input,fxyz_cn1_lc3,out
/input,fxyz_cn2_lc3,out
/input,fxyz_cn3_lc3,out

/com Distribute TF Loads according to node location
xl1=1.8194        $yl1=0.79397     $zl1=-1.2319     $fl1=  3700.00000  !b outboard bottom
xl2=1.847          $yl2=0.77363     $zl2=1.2319      $fl2= -3700.00000  !b outboard top
xl3=0.92238      $yl3=0.25725     $zl3=1.3767      $fl3= -17000.0000  !A inboard top left
xl4=0.96987      $yl4=0.19198     $zl4=-1.3767     $fl4=  17000.0000  !A inboard bottom left  
xl5=1.9958        $yl5=0.13282     $zl5=-1.2319     $fl5=  2100.00000  !a outboard bottom left 
xl6=1.9216        $yl6=0.55548     $zl6=-1.2319     $fl6=  4800.00000  !a outboard bottom right 
xl7=1.9958        $yl7=0.1328      $zl7=1.2319      $fl7= -2100.00000  !A outboard top left  
xl8=1.9216        $yl8=0.55548     $zl8=1.2319      $fl8= -4800.00000  !A outboard top right  
xl9=0.6203        $yl9=0.75607     $zl9=1.1986      $fl9= -17000.0000  !c inboard top 
xl10=1.4484      $yl10=1.3887     $zl10=-1.1986    $fl10= 7300.00000  !C outboard bottom left 
xl11=1.1155      $yl11=1.6664     $zl11=-1.1986    $fl11= 3900.00000  !C outboard bottom right  
xl12=1.4482      $yl12=1.3886     $zl12=1.1986     $fl12=-7300.00000  !C outboard top left
xl13=1.1691      $yl13=1.6298     $zl13=1.1986     $fl13=-3900.00000  !C outboard top right
xl14=0.61696     $yl14=0.75987   $zl14=-1.1986    $fl14= 17000.0000  !C inboard bottom
            
allsel
csys
*do,i,1,14
fn=node(xl%i%,yl%i%,zl%i%)
f,fn,FZ,fl%i%      
*enddo

allsel

!esel,u,num,,314079
!esel,u,num,,314096
!esel,u,num,,314113
!esel,u,num,,314130
!esel,u,num,,314147
!esel,u,num,,314164

nsubst,10,20,2
outres,all,1
solve

!/solu
!ANTYPE,,REST, , ,0  
!DOFSEL,S,FX,FY,FZ   
!FSCALE,1e-4,
!DOFSEL,ALL  
!solve   

FINI
/post1
set,last
/com
/com Plot bolt loads
/com
esel,s,type,,101,104
nsle
/com Axial Member Forces
etab,fxi,smisc,1
etab,fxj,smisc,7
sadd,fxikip,fxi,,0.2248/1000
sadd,fxjkip,fxj,,0.2248/1000
/com Transverse Member Forces
etab,fy_i,smisc,2
etab,fy_j,smisc,8
etab,fz_i,smisc,3
etab,fz_j,smisc,9
smult,fyi2,fy_i,fy_i
smult,fzi2,fz_i,fz_i
sadd,fyzi2,fyi2,fzi2
sexp,fyzi,fyzi2,,0.5
smult,fyj2,fy_j,fy_j
smult,fzj2,fz_j,fz_j
sadd,fyzj2,fyj2,fzj2
sexp,fyzj,fyzj2,,0.5
sadd,fyzikip,fyzi,,0.2248/1000
sadd,fyzjkip,fyzj,,0.2248/1000
/com Print Member Forces
pret,fxikip,fyzikip
/com Plot Member Forces
!PLLS,FXI,FXJ,1,0  
PLLS,FXIkip,FXJkip,1,0
/capture,tensile,jpg  
!PLLS,FYZI,FYZJ,1,0  
PLLS,FYZIkip,FYZJkip,1,0 
/capture,shear, 

/com Contact Interface
!esel,s,type,,49     ! A-B
!esel,r,real,,49     ! A-B
!esel,s,type,,55     ! BB-C
!esel,r,real,,55     ! BB-C
!esel,s,type,,66     ! A-A
!esel,r,real,,66     ! A-A
!esel,s,type,,64     ! C-C
!esel,r,real,,64     ! C-C
!nsle
!etab,cpres,smisc,13
!plet,cpres,avg
!etab,aasri,nmisc,50
!etab,aasrj,nmisc,51
!etab,aasrk,nmisc,52
!etab,aasrl,nmisc,53
!sadd,aasr1,aasri,aasrj
!sadd,aasr2,aasrk,aasrl
!sadd,aasr,aasr1,aasr2,0.25,0.25
!etab,aassi,nmisc,54
!etab,aassj,nmisc,55
!etab,aassk,nmisc,56
!etab,aassl,nmisc,57
!sadd,aass1,aassi,aassj
!sadd,aass2,aassk,aassl
!sadd,aass,aass1,aass2,0.25,0.25
!smult,aasr2,aasr,aasr
!smult,aass2,aass,aass
!sadd,aa2,aasr2,aass2
!sexp,aa,aa2,,0.5
!sadd,slip_mil,aa,,1000/0.0254
!pret,aasri,aasrj,aasrk,aasrl,aasr
!pret,aassi,aassj,aassk,aassk,aass
!pret,aasr,aass,aa
!plet,slip_mil,avg

*if,k_bolts,eq,1,then
esel,s,type,,66
nsle
*elseif,k_bolts,eq,2
esel,s,type,,49
nsle
*elseif,k_bolts,eq,3
esel,s,type,,55
nsle
*elseif,k_bolts,eq,4
esel,s,type,,64
nsle
*endif

PLESOL, CONT,PRES, 0,1.0
PLESOL, CONT,SLIDE, 2,1.0   
PLESOL, CONT,STAT, 0,1.0

fini
/exit,all
/eof

/com Calc Normal Stress 
esel,s,mat,,16
nsle
csys,1
rsys,1
esel,s,mat,,16
nsle
etab,sth,s,y
etab,v,volu
taa=0.00635   ! thickness of shim
smult,sthv_a,sth,v,1/taa,0.2248
ssum
*get,fth_lb,ssum,,item,sthv_a

/com
/com Integrate Shears over phi-bands at each interface
/com
:100
local,2003,3,,,,,,,2
*dim,shears,,100,8,4
/com
/com AA: Shim:mat16, CONT:type66/real65,
/com
phiaamn=145   ! angular start of inboard leg
phiaamx=215   ! angular   end of inboard leg
dphiaa=2.5    ! angular extent of integration
taa=0.00635   ! thickness of shim
esel,s,mat,,16
nsle
rsys,1
etab,sxy,s,xy
etab,syz,s,yz
etab,sy,s,y
etab,v,volu
smult,fxyv,sxy,v,,1/taa
smult,fyzv,syz,v,,1/taa
smult,fyv,sy,v,,1/taa
csys,2003
nsel,r,loc,z,phiaamn,phiaamx
esln,,1
esel,r,mat,,16
cm,e_ibaa,elem
esel,none
cm,e_done,elem
*do,j,1,nint((phiaamx-phiaamn)/dphiaa)
cmsel,s,e_ibaa
nsle
nsel,r,loc,z,phiaamn+(j-1.5)*dphiaa,phiaamn+(j+0)*dphiaa
shears(j,1,1)=(phiaamn+(j-1)*dphiaa+phiaamn+(j+0)*dphiaa)/2
esln,,1
esel,r,mat,,16
cmsel,u,e_done
ssum
*get,shears(j,3,1),ssum,,item,fxyv                    ! Radial Shear Force, N
*get,shears(j,4,1),ssum,,item,fyzv                    ! Vertical Shear Force, N
*get,shears(j,8,1),ssum,,item,fyv                     ! Normal Force, N

*get,vj,ssum,,item,v
shears(j,2,1)=vj/taa                                  ! Area of Patch, m^2
shears(j,5,1)=shears(j,3,1)/shears(j,2,1)/1e6         ! Stress, MPa
shears(j,6,1)=shears(j,4,1)/shears(j,2,1)/1e6         ! Stress, MPa
shears(j,7,1)=sqrt(shears(j,5,1)**2+shears(j,6,1)**2) ! Stress, MPa
cmsel,a,e_done
cm,e_done,elem
*enddo


/output,shearsaa2%rn%a,out
*vwrite,
('      Angle(deg)     Area(m**2)       Fxy(N)         Fyz(N)        Txy(MPa)       Tyz(MPa)       Tsrs(MPa)      Fy(N)')
*vwri,shears(1,1,1),shears(1,2,1),shears(1,3,1),shears(1,4,1),shears(1,5,1),shears(1,6,1),shears(1,7,1),shears(1,8,1)
(1x,1p8e15.4)
/output
/eof

:200
/com
/com AB: Shim:mat13, CONT:type49/real50,
/com
phiabmn=145   ! angular start of inboard leg
phiabmx=215   ! angular   end of inboard leg
dphiab=2.5    ! angular extent of integration
tab=0.0127    ! thickness of shim
esel,s,mat,,13
nsle
rsys,1
etab,sxy,s,xy
etab,syz,s,yz
etab,v,volu
smult,fxyv,sxy,v,,1/tab
smult,fyzv,syz,v,,1/tab
csys,2003
nsel,r,loc,z,phiabmn,phiabmx
esln,,1
esel,r,mat,,13
cm,e_ibab,elem
esel,none
cm,e_done,elem
*do,j,1,nint((phiabmx-phiabmn)/dphiab)
cmsel,s,e_ibab
nsle
nsel,r,loc,z,phiabmn+(j-1.5)*dphiab,phiabmn+(j+0)*dphiab
shears(j,1,2)=(phiabmn+(j-1)*dphiab+phiabmn+(j+0)*dphiab)/2
esln,,1
esel,r,mat,,13
cmsel,u,e_done
ssum
*get,shears(j,3,2),ssum,,item,fxyv                    ! Radial Shear Force, N
*get,shears(j,4,2),ssum,,item,fyzv                    ! Vertical Shear Force, N
*get,vj,ssum,,item,v
shears(j,2,2)=vj/tab                                  ! Area of Patch, m^2
shears(j,5,2)=shears(j,3,2)/shears(j,2,2)/1e6         ! Stress, MPa
shears(j,6,2)=shears(j,4,2)/shears(j,2,2)/1e6         ! Stress, MPa
shears(j,7,2)=sqrt(shears(j,5,2)**2+shears(j,6,2)**2) ! Stress, MPa
cmsel,a,e_done
cm,e_done,elem
*enddo


/output,shearsab2%rn%,out
*vwrite,
('      Angle(deg)     Area(m**2)       Fxy(N)         Fyz(N)        Txy(MPa)       Tyz(MPa)       Tsrs(MPa)')
*vwri,shears(1,1,2),shears(1,2,2),shears(1,3,2),shears(1,4,2),shears(1,5,2),shears(1,6,2),shears(1,7,2)
(1x,1p7e15.4)
/output

:300
/com
/com BC: Shim:mat14, CONT:type55/real56
/com
phibcmn=145   ! angular start of inboard leg
phibcmx=205   ! angular   end of inboard leg
dphibc=2.5    ! angular extent of integration
tbc=0.0127    ! thickness of shim
esel,s,mat,,14
nsle
rsys,1
etab,sxy,s,xy
etab,syz,s,yz
etab,v,volu
smult,fxyv,sxy,v,,1/tbc
smult,fyzv,syz,v,,1/tbc
csys,2003
nsel,r,loc,z,phibcmn,phibcmx
esln,,1
esel,r,mat,,14
cm,e_ibbc,elem
esel,none
cm,e_done,elem
*do,j,1,nint((phibcmx-phibcmn)/dphibc)
cmsel,s,e_ibbc
nsle
nsel,r,loc,z,phibcmn+(j-1.5)*dphibc,phibcmn+(j+0)*dphibc
shears(j,1,3)=(phibcmn+(j-1)*dphibc+phibcmn+(j+0)*dphibc)/2
esln,,1
esel,r,mat,,14
cmsel,u,e_done
ssum
*get,shears(j,3,3),ssum,,item,fxyv                    ! Radial Shear Force, N
*get,shears(j,4,3),ssum,,item,fyzv                    ! Vertical Shear Force, N
*get,vj,ssum,,item,v
shears(j,2,3)=vj/tbc                                  ! Area of Patch, m^2
shears(j,5,3)=shears(j,3,3)/shears(j,2,3)/1e6         ! Stress, MPa
shears(j,6,3)=shears(j,4,3)/shears(j,2,3)/1e6         ! Stress, MPa
shears(j,7,3)=sqrt(shears(j,5,3)**2+shears(j,6,3)**2) ! Stress, MPa
cmsel,a,e_done
cm,e_done,elem
*enddo


/output,shearsbc2%rn%,out
*vwrite,
('      Angle(deg)     Area(m**2)       Fxy(N)         Fyz(N)        Txy(MPa)       Tyz(MPa)       Tsrs(MPa)')
*vwri,shears(1,1,3),shears(1,2,3),shears(1,3,3),shears(1,4,3),shears(1,5,3),shears(1,6,3),shears(1,7,3)
(1x,1p7e15.4)
/output


:400
/com
/com CC: Shim:mat15, CONT:type64/real66
/com
phiccmn=140   ! angular start of inboard leg
phiccmx=220   ! angular   end of inboard leg
dphicc=2.5    ! angular extent of integration
tcc=0.00635   ! thickness of shim
esel,s,mat,,15
nsle
rsys,1
etab,sxy,s,xy
etab,syz,s,yz
etab,v,volu
smult,fxyv,sxy,v,,1/tcc
smult,fyzv,syz,v,,1/tcc
csys,2003
nsel,r,loc,z,phiccmn,phiccmx
esln,,1
esel,r,mat,,15
cm,e_ibcc,elem
esel,none
cm,e_done,elem
*do,j,1,nint((phiccmx-phiccmn)/dphicc)
cmsel,s,e_ibcc
nsle
nsel,r,loc,z,phiccmn+(j-1.5)*dphicc,phiccmn+(j+0)*dphicc
shears(j,1,4)=(phiccmn+(j-1)*dphicc+phiccmn+(j+0)*dphicc)/2
esln,,1
esel,r,mat,,15
cmsel,u,e_done
ssum
*get,shears(j,3,4),ssum,,item,fxyv                    ! Radial Shear Force, N
*get,shears(j,4,4),ssum,,item,fyzv                    ! Vertical Shear Force, N
*get,vj,ssum,,item,v
shears(j,2,4)=vj/tcc                                  ! Area of Patch, m^2
shears(j,5,4)=shears(j,3,4)/shears(j,2,4)/1e6         ! Stress, MPa
shears(j,6,4)=shears(j,4,4)/shears(j,2,4)/1e6         ! Stress, MPa
shears(j,7,4)=sqrt(shears(j,5,4)**2+shears(j,6,4)**2) ! Stress, MPa
cmsel,a,e_done
cm,e_done,elem
*enddo


/output,shearscc2%rn%,out
*vwrite,
('      Angle(deg)     Area(m**2)       Fxy(N)         Fyz(N)        Txy(MPa)       Tyz(MPa)       Tsrs(MPa)')
*vwri,shears(1,1,4),shears(1,2,4),shears(1,3,4),shears(1,4,4),shears(1,5,4),shears(1,6,4),shears(1,7,4)
(1x,1p7e15.4)
/output


fini
/EXIT,ALL   
/EOF


!!!!post processing
!!! produces stress images along the path of the winding normal to the tee.


/prep7

!!!! plots stresses near clamps on coil c
cl_num = 35    !A coil =     B Coil = 9 through 114 C coil = 15 through 920
*do,i,889,890    ! A coil = 829-876, B coil = 68-116,  C Coil = 879-926

csys,i
wpcsys
wprota,,90,90  !exclude for coil 2, ,,-90,90 for coil A
/cplane,1
/type,,1
/view,,wp
/dist,,.25   !.4 for far view, .1 for near
/replot
cswpla,5000
nsel,s,loc,z,-.0254,.0254 	!select one inch deep
nsel,r,loc,x,-.5,.5 		!select out IN X (.5 FOR FAR, .1 FOR NEAR)
nsel,r,loc,y,-.5,.5		!select out IN Y (.5 FOR FAR, .1 FOR NEAR)
esln
/post1
jpeg,qual,100
esel,s,type,,7
!esel,u,name,,badclamp
!esel,a,name,,clamp3
rsys,5000
/plopts,info,auto
/plopts,leg2,off
/plopts,date,off
/plopts,wp,on

!!!plots three x, y, z stresses

!/Title,Von_Mises Stress near clamp %cl_num% on coil A
!/dscale,,off
!plnsol,s,eqv,0,1
!/image,save,Clamp_%cl_num%_Von_Mises_stress,bmp
!/ui,copy,save,jpeg,full,color,norm,portrait,no,100
!/rename,non_linear000,jpg,,Hole_%Cl_num%_Stres_Intensity,jpg

/Title,Y Stress near clamp %cl_num% on coil A
/dscale,,off
/contour,,18,-.14e9,,.14e9
plnsol,s,Y,0,1
/ui,copy,save,jpeg,full,color,norm,portrait,no,100
/rename,non_linear000,jpg,,Y_stress_Hole_%Cl_num%,jpg

/Title,Z Stress near clamp %cl_num% on coil A
/dscale,,off
/contour,,18,-.14e9,,.14e9
plnsol,s,Z,0,1
/ui,copy,save,jpeg,full,color,norm,portrait,no,test
/rename,non_linear000,jpg,,Z_stress_Hole_%Cl_num%,jpg

/Title,x Stress near clamp %cl_num% on coil A
/dscale,,off
/contour,,18,-.14e9,,.14e9
plnsol,s,x,0,1
/ui,copy,save,jpeg,full,color,norm,portrait,yes,100
/rename,non_linear000,jpg,,X_stress_Hole_%Cl_num%,jpg

!!!! plot deformations

!/Title,X deflection near clamp %cl_num%
!/dscale,,5
!plnsol,u,X,2
!/image,save,X_deflection_near_clamp_%cl_num%,bmp

!/Title,Y deflection near clamp %cl_num%
!/dscale,,5
!plnsol,u,Y,2
!/image,save,Y_deflection_near_clamp_%cl_num%,bmp

!/Title,Z deflection near clamp %cl_num%
!/dscale,,5
!plnsol,u,Z,2
!/image,save,Z_deflection_near_clamp_%cl_num%,bmp

cl_num = cl_num + 2
*enddo



/com
/com
/com Select all nodes on A-side holes
nsel,none
*do,j,1,26
cmsel,a,holea%j%
*enddo

/com Select all nodes on A-side holes
nsel,none
*do,j,1,26
cmsel,a,holeb%j%
*enddo



/eshap,1,1
plet,si_ij,avg