然后现在是建模的原始文件：
def in_def
xx = in(msg+'('+'default:'+string(default)+'):')
if type(xx) = 3
in_def = default
else
in_def = xx
endif
end
;
def beam_canshu
default = 6.0
msg=' yuan gong banjing (r) '
r = in_def
a13=r
b13=0
a1=-a13
b1=0
a12=r*cos(pi/12)
b12=r*sin(pi/12)
a2=-a12
b2=b12
a11=r*cos(pi/6)
b11=r*sin(pi/6)
a3=-a11
b3=b11
a10=r*cos(pi/4)
b10=r*sin(pi/4)
a4=-a10
b4=b10
a9=r*cos(pi/3)
b9=r*sin(pi/3)
a5=-a9
b5=b9
a8=r*cos((pi*5)/12)
b8=r*sin((pi*5)/12)
a6=-a8
b6=b8
a7=0
b7=r
;
default = -5e4
msg=' cross force thickness1 (p0) '
p0 = in_def
;
default = -1e5
msg=' cross force xielv (l0) '
l0 = in_def
;
default = 0
msg=' vertical force xielv (c0) '
c0 = in_def
;
default = 5e04
msg=' horizontal crustal stress xiaojianju (m0) '
m0 = in_def
;
default = 1e06
msg=' horizontal crustal stress yuanlice (n0) '
n0 = in_def
;
f1=(p0+l0*r)*r*sin(pi/12)
f2=((p0+l0*r)*r*sin(pi/6))/2+(l0*sin(pi/12)*2)/3+(l0*(sin(pi/6)-sin(pi/12)))/3
f3=((p0+l0*r)*r*(sin(pi/4)-sin(pi/12)))/2+(l0*(sin(pi/6)-sin(pi/12))*2)/3+(l0*(sin(pi/4)-sin(pi/6)))/3
f4=((p0+l0*r)*r*(sin(pi/3)-sin(pi/6)))/2+(l0*(sin(pi/4)-sin(pi/6))*2)/3+(l0*(sin(pi/3)-sin(pi/4)))/3
f5=((p0+l0*r)*r*(sin((pi*5)/12)-sin(pi/4)))/2+(l0*(sin(pi/3)-sin(pi/4))*2)/3+(l0*(sin((pi*5)/12)-sin(pi/3)))/3
f6=((p0+l0*r)*r*(1-sin(pi/3)))/2+(l0*(sin((pi*5)/12)-sin(pi/3))*2)/3+(l0*(1-sin((pi*5)/12)))/3
f7=((p0+l0*r)*r*sin(pi/6))/2-(l0*sin(pi/12)*2)/3-(l0*(sin(pi/6)-sin(pi/12)))/3
f8=((p0+l0*r)*r*(sin(pi/4)-sin(pi/12)))/2-(l0*(sin(pi/6)-sin(pi/12))*2)/3-(l0*(sin(pi/4)-sin(pi/6)))/3
f9=((p0+l0*r)*r*(sin(pi/3)-sin(pi/6)))/2-(l0*(sin(pi/4)-sin(pi/6))*2)/3-(l0*(sin(pi/3)-sin(pi/4)))/3
f10=((p0+l0*r)*r*(sin((pi*5)/12)-sin(pi/4)))/2-(l0*(sin(pi/3)-sin(pi/4))*2)/3-(l0*(sin((pi*5)/12)-sin(pi/3)))/3
f11=((p0+l0*r)*r*(1-sin(pi/3)))/2-(l0*(sin((pi*5)/12)-sin(pi/3))*2)/3-(l0*(1-sin((pi*5)/12)))/3
;
c1=(c0*(1-cos(pi/12))*2+c0*(cos(pi/6)-cos(pi/12)))/3+(n0)*r(1-cos(pi/6))
c2=(c0*(cos(pi/6)-cos(pi/12))*2+c0*(cos(pi/4)-cos(pi/6)))/3+(n0)*r*(cos(pi/4)-cos(pi/12))
c3=(c0*(cos(pi/4)-cos(pi/6))*2+c0*(cos(pi/3)-cos(pi/4)))/3+(n0)*r*(cos(pi/3)-cos(pi/6))
c4=(c0*(cos(pi/3)-cos(pi/4))*2+c0*(cos((pi*5)/12)-cos(pi/3)))/3+(n0)*r*(cos((pi*5)/12)-cos(pi/4))
c5=(c0*(cos((pi*5)/12)-cos(pi/3))*2+c0*(1-cos((pi*5)/12)))/3+(n0)*r*cos(pi/3)
c6=-((c0*(cos((pi*5)/12)-cos(pi/3))*2+c0*(1-cos((pi*5)/12)))/3+(m0)*r*cos(pi/3))
c7=-((c0*(cos(pi/3)-cos(pi/4))*2+c0*(cos((pi*5)/12)-cos(pi/3)))/3+(m0)*r*(cos((pi*5)/12)-cos(pi/4)))
c8=-((c0*(cos(pi/4)-cos(pi/6))*2+c0*(cos(pi/3)-cos(pi/4)))/3+(m0)*r*(cos(pi/3)-cos(pi/6)))
c9=-((c0*(cos(pi/6)-cos(pi/12))*2+c0*(cos(pi/4)-cos(pi/6)))/3+(n0)*r*(cos(pi/4)-cos(pi/12)))
c10=-((c0*(1-cos(pi/12))*2+c0*(cos(pi/6)-cos(pi/12)))/3+(n0)*r(1-cos(pi/6)))
cc=((n0-m0)*(1-cos(pi/12)))/2
end
beam_canshu
sel beam id=1 begin=( a1,b1,0 ) end=( a2,b2,0 ) nesg=1
sel beam id=1 begin=( a2,b2,0 ) end=( a3,b3,0 ) nesg=1
sel beam id=1 begin=( a3,b3,0 ) end=( a4,b4,0 ) nesg=1
sel beam id=1 begin=( a4,b4,0 ) end=( a5,b5,0 ) nesg=1
sel beam id=1 begin=( a5,b5,0 ) end=( a6,b6,0 ) nesg=1
sel beam id=1 begin=( a6,b6,0 ) end=( a7,b7,0 ) nesg=1
sel beam id=1 begin=( a7,b7,0 ) end=( a8,b8,0 ) nesg=1
sel beam id=1 begin=( a8,b8,0 ) end=( a9,b9,0 ) nesg=1
sel beam id=1 begin=( a9,b9,0 ) end=( a10,b10,0 ) nesg=1
sel beam id=1 begin=( a10,b10,0 ) end=( a11,b11,0 ) nesg=1
sel beam id=1 begin=( a11,b11,0 ) end=( a12,b12,0 ) nesg=1
sel beam id=1 begin=( a12,b12,0 ) end=( a13,b13,0 ) nesg=1
plot set rot 90 0 0
sel beam id=1 prop emod=52.4e9 nu=0.30 xcarea=3.55e-3 xcj=400e-6 xciy=333e-6 xciz=67e-6
sel node fix z xr yr
sel node fix x y zr range id=1
sel node fix x y zr range id=13
sel node apply force=(cc,f1,0.0) range id=7
sel node apply force=(c1,f2,0.0) range id=8
sel node apply force=(c2,f3,0.0) range id=9
sel node apply force=(c3,f4,0.0) range id=10
sel node apply force=(c4,f5,0.0) range id=11
sel node apply force=(c5,f6,0.0) range id=12
sel node apply force=(c6,f11,0.0) range id=2
sel node apply force=(c7,f10,0.0) range id=3
sel node apply force=(c8,f9,0.0) range id=4
sel node apply force=(c9,f8,0.0) range id=5
sel node apply force=(c10,f7,0.0) range id=6
his sel node ydisp id=7
his sel beam force fx end1 cid=2
his sel beam force fx end1 cid=4
his sel beam force fx end1 cid=6
his sel beam force fx end1 cid=7
his sel beam force fx end1 cid=9
his sel beam force fx end1 cid=11
plot create geom_and_disp
plot set title text
Geometry of the modeled system
plot set cap size 25
plot set rot 90 0 0
plot add sel geom black black
plot add axes black
plot add sel disp blue
plot create shear_force
plot set title text
Shear force distribution in the beam
plot set cap size 25
plot set rot 90 0 0
plot add sel geom black black node=off
plot add sel beam force fy scale=0.06
plot create axis_force
plot set title text
axis force in the beam
plot set cap size 25
plot set rot 90 0 0
plot add sel geom black black node=off
plot add sel beam force fx
plot create moment
plot set title text
Bending moment distribution in the beam
plot set cap size 25
plot set rot 90 0 0
plot add sel geom black black node=off
plot add sel beam moment mz scale=0.06
print sel beam id
print sel beam local
print sel beam position
solve ratio=1e-7
print sel beam force
