*SET,ALPH,0.5
*SET,TEMP,1
a=100
c1=0.4988
c2=1-c1
r1=sqrt(c1*a*a/3.1415926*4)
ET,1,PLANE42
KEYOPT,1,3,2
MP,EX,1,83.3
MP,PRXY,1,0.22
MPTEMP,,,,,,,,
MPTEMP,1,0
UIMP,1,REFT,,,
MPDATA,ALPX,1,,ALPH
MPDATA,ALPY,1,,-ALPH
MPDATA,ALPZ,1,,0
MP,EX,2,3.33
MP,PRXY,2,0.35
MPTEMP,,,,,,,,
MPTEMP,1,0
UIMP,1,REFT,,,
MPDATA,ALPX,2,,ALPH
MPDATA,ALPY,2,,-ALPH
MPDATA,ALPZ,2,,0
RECTNG,0,a,0,a,
PCIRC,r1, ,0,90,
AOVLAP,all
wpro,-45.000000,,
wpro,,,-90.000000
asbw,4
WPCSYS,-1,0
WPROTA,-45
CSWPLA,11,0,1,1,
CSYS,11
lsel,s,,,2,4
lsel,a,,,6
LESIZE,ALL, , ,11, ,1, , ,1,
lsel,s,,,10,11
lsel,a,,,1
LESIZE,ALL, , ,6, ,1, , ,1,
lsel,s,,,8,9
LESIZE,ALL, , ,22, ,1, , ,1,
allsel,
TYPE,1
MAT,1
ESYS,11
MSHAPE,0,2D
MSHKEY,0
amesh,3
TYPE,1
MAT,2
ESYS,11
MSHAPE,0,2D
MSHKEY,1
amesh,1,2
WPROTA,45
CSYS,0
ARSYM,X,ALL, , , ,0,0
ARSYM,Y,ALL, , , ,0,0
NUMMRG,ALL
NSEL,S,LOC,Y,-a
NSEL,a,LOC,Y,a
d,all,,,,,,uX
NSEL,S,LOC,x,a
NSEL,a,LOC,x,-a
d,all,,,,,,uy
NSEL,ALL
TUNIF,TEMP,
ALLSEL,ALL
/SOLU
SOLVE
/POST1
rsys,0
SET,LAST
KA=1/(TEMP*ALPH*2)
ESEL,ALL
*get,MM,elem,,num,min
*get,coun,elem,,count
AREA=0
C=0
*do,I,1,coun
*get,a1,elem,MM,area
*get,MM,elem,MM,nxth
area=area+a1
*enddo
*DIM,S,,4
*DIM,T,,4
*DIM,N,,4
*DIM,X,,4
*DIM,Y,,4
*DIM,U,,4
*DIM,V,,4
*DIM,NS,,4
*DIM,NT,,4
*DIM,NIX,,4
*DIM,NIY,,4
S(1)=0.577350269189626
S(2)=S(1)
S(3)=-S(1)
S(4)=-S(1)
T(1)=S(1)
T(2)=-S(1)
T(3)=S(1)
T(4)=-S(1)
E1212H=0
*do,sel,1,2
ESEL,,MAT,,sel
*GET,M,ELEM,,NUM,MIN
*GET,COUN,ELEM,,COUNT
*GET,EPSI,EX,sel
*GET,NU,NUXY,sel
*DO,I,1,COUN !对每种材料的单元数目进行循环
PX=0
PY=0
N(1)=NELEM(M,1) !***********单元节点编号
N(2)=NELEM(M,2)
N(3)=NELEM(M,3)
N(4)=NELEM(M,4) !***********单元节点编号
X(1)=NX(N(1)) !***********各节点坐标
X(2)=NX(N(2))
X(3)=NX(N(3))
X(4)=NX(N(4))
Y(1)=NY(N(1))
Y(2)=NY(N(2))
Y(3)=NY(N(3))
Y(4)=NY(N(4)) !***********各节点坐标
U(1)=UX(N(1)) !***********各节点横向,竖向位移
U(2)=UX(N(2))
U(3)=UX(N(3))
U(4)=UX(N(4))
V(1)=UY(N(1))
V(2)=UY(N(2))
V(3)=UY(N(3))
V(4)=UY(N(4)) !***********各节点横向,竖向位移
*DO,J,1,4 !
NS(1)=-0.25*(1-T(J)) !***********计算插值点的形函数
NS(2)=0.25*(1-T(J))
NS(3)=0.25*(1+T(J))
NS(4)=-0.25*(1+T(J))
NT(1)=-0.25*(1-S(J))
NT(2)=-0.25*(1+S(J))
NT(3)=0.25*(1+S(J))
NT(4)=0.25*(1-S(J)) !***********计算插值点的形函数
XS=0
XT=0
YS=0
YT=0
*DO,K,1,4 !***********计算XS,XT,YS,YT,坐标变换式
XS=NS(K)*X(K)+XS
XT=NT(K)*X(K)+XT
YS=NS(K)*Y(K)+YS
YT=NT(K)*Y(K)+YT
*ENDDO !***********计算XS,XT,YS,YT,坐标变换式
*DO,K,1,4 !***********计算NIX,NIY,应变矩阵元素
NIX(K)=YT*NS(K)-YS*NT(K)
NIY(K)=-XT*NS(K)+XS*NT(K)
*ENDDO !***********计算NIX,NIY,应变矩阵元素
EPX=0
EPY=0
*DO,K,1,4 !***********计算插值点应变
EPX=NIX(K)*V(K)+EPX
EPY=NIY(K)*U(K)+EPY
*ENDDO !***********计算插值点应变
PX=PX+EPX
PY=PY+EPY
*ENDDO !***********计算单元应变PX,PY
*GET,DA,ELEM,M,AREA !提取单元的面积
K=EPSI/(1+NU)/2/AREA
E1212H=K*(DA-PY*KA-PX*KA
*GET,M,ELEM,M,NXTH
*ENDDO
*enddo
allsel,
finish
G12=E1212H