finish
/clear,nostart
/units,si

	!!!material properties
density=7850		!density
young=2.1e11		!Young modulus
poisson=0.3		!Poisson's ratio
CTEX=11.1e-6		!thermal expansion coefficient
g=9.81			!gravity acceleration

	!!!geometry
a=40e-3			!width (X-axis)
l=40e-3			!length (Z-axis)
h=60e-3			!heigth (Y-axis)
w=20e-3			!wall thickness

	!!!queries
*ask,case,Choose case (A/B/C) enclosed in single quotes,'a'
*if,case,eq,'A',or,case,eq,'a',then
	*ask,f,Enter Force (N),160e3
	*ask,dt,Enter tempeature difference (K),100
*elseif,case,eq,'B',or,case,eq,'b',then
	*ask,dt,Enter tempeature difference (K),100
*else
	*ask,f,Enter Force (N),160e3
	*ask,mu,Enter friction coefficient,0.12
*endif
*ask,na,Enter number of divisions for width,10
*ask,nl,Enter number of divisions for length,10
*ask,nh,Enter number of divisions for heigth,10
*ask,nw,Enter number of divisions for wall thickness,3

/prep7
	!!!element types
et,1,solid186
et,2,targe170
et,3,conta174
r,1
rmore,0,0

	!!!material properties
mp,dens,1,density
mp,ex,1,young
mp,nuxy,1,poisson
uimp,1,reft,,,0
mp,ctex,1,ctex
*if,mu,ne,0,then
	mp,mu,1,mu
*endif

	!!!slot geometry
block,-a/2-w,a/2+w,-w,0,-l/2-w,l/2+w
block,-a/2-w,-a/2,-w,h,-l/2-w,l/2+w
block,a/2,a/2+w,-w,h,-l/2-w,l/2+w
block,-a/2-w,a/2+w,-w,h,-l/2-w,-l/2
block,-a/2-w,a/2+w,-w,h,l/2,l/2+w
vovlap,all
cm,slot,volu		!slot geometry extracted to volume component

	!!!parallelipiped
block,-a/2,a/2,0,h,-l/2,l/2
vsel,u,volu,,slot
*get,v1,volu,,num,min	!specimen volume number stored in v1 variable
allsel

	!!!meshing
mshape,0,3d
mshkey,1
lsel,s,loc,x,0
lesize,all,,,na
lsel,s,loc,z,0
lesize,all,,,nl
lsel,s,loc,y,h/2
lesize,all,,,nh
allsel
lesize,all,,,nw
vmesh,all
vsel,s,volu,,v1
aslv,s
asel,u,loc,y,h
nsla,s,1
type,3
esurf
asel,s,loc,x,-a/2,a/2
asel,r,loc,z,-l/2,l/2
asel,r,loc,y,0,h
vsel,inve
aslv,r
nsla,s,1
type,2
esurf
allsel

	!!!boundary conditions
acel,,g
tunif,0
*if,case,eq,'A',or,case,eq,'a',then
	vsel,s,volu,,v1
	aslv,s
	asel,r,loc,y,h
	nsla,s,1
	sf,all,pres,f/(a*l)
	vsel,s,volu,,v1
	allsel,below,volu
	bfe,all,temp,1,dt
	vsel,s,volu,,slot
	allsel,below,volu
	d,all,all
*elseif,case,eq,'B',or,case,eq,'b',then
	vsel,s,volu,,v1
	allsel,below,volu
	bfe,all,temp,1,dt
	vsel,s,volu,,slot
	allsel,below,volu
	d,all,all
*else
	vsel,s,volu,,v1
	aslv,s
	asel,r,loc,y,h
	nsla,s,1
	sf,all,pres,f/(a*l)
*endif
allsel
asel,s,loc,y,-w
asel,a,loc,x,-a/2-w
asel,a,loc,x,a/2+w
asel,a,loc,z,-l/2-w
asel,a,loc,z,l/2+w
nsla,s,1
d,all,all

allsel

/solu

solve

/post1



/view,1,1,1,1

esel,s,type,,3
etable,stot,cont,stot  
esel,s,cent,x,a/2
esel,a,cent,x,-a/2
esel,r,type,,3
ssum
*get,stotx,ssum,,item,stot
esel,s,cent,z,l/2
esel,a,cent,z,-l/2
esel,r,type,,3
ssum
*get,stotz,ssum,,item,stot
esel,s,cent,y,0
esel,r,type,,3
ssum
*get,stoty,ssum,,item,stot

set,last
esel,s,cent,x,0,a/2+w
esel,r,cent,z,0,l/2+w
esel,inve

*if,case,eq,'A',or,case,eq,'a',then
/title, Case A with temperature load: F = %f% N, dT = %dt%
*elseif,case,eq,'B',or,case,eq,'b',then
/title, Case B: dT = %dt%
*else
/title, Case C: F = %f% N, mu = %mu%
*endif

plnsol,u,y,2
sx=stotx/nh/nl/2
sy=stoty/nl/na
sz=stotz/na/nh/2
/tlabel,-1,-.6,Total stresses on contact surfaces in X-direction %sx% Pa
/tlabel,-1,-.65,Total stresses on contact surfaces in Y-direction %sy% Pa
/tlabel,-1,-.7,Total stresses on contact surfaces in Z-direction %sz% Pa