- 积分
- 0
- 注册时间
- 2011-12-21
- 仿真币
-
- 最后登录
- 1970-1-1
|
就是英文手册后面的一个例子。主要是里面定义的fish语言看不懂。
;---------------------------------------------------------------------
;triax.dat compression test of cylindrical sample using
; ubiquitous joint model
;---------------------------------------------------------------------
def test0
tim0 = clock
end
def tim
tim = (clock-tim0)/100.0
end
gen zone brick p0 0 0 0 p1 15 0 0 p2 0 5 0 p3 0 0 15 size 15 5 15
;
def p_cons
eps = 0.001
rad = 1.
c_area = 3.
betaj = 1
if beta = 0. then
betaj = 0
end_if
if beta = 90. then
betaj = 0
end_if
end
;
def trisol
loop k (0,18)
kk = k
beta = k * 15.
p_cons
command
model null
model ubiquitous
prop bul 1.e8 shea 7.e7 cohesion 2.e3
prop friction 40. dilation 0. tension 2400.
prop jdip beta jdd 0. jcoh 1.e3 jfric 30. jdil 0. jten 2000.
fix z range z -.1 .1
fix z range z 14.9 15.1
ini zvel 1.e-7 range z -.1 .1
ini zvel -1.e-7 range z 14.9 15.1
cyc 4500
print beta
print sigmav
print anasol
print p_err
print tim
end_command
end_loop
end
;
;---------------------------------------------------------------------
; sigmav : numerical value of the compressive strength
; p_sol : analytical prediction for the compressive strength
; p_err : relative difference between sigmav and p_sol (in %)
;---------------------------------------------------------------------
def anasol
pnt = zone_head
mc = z_prop(pnt,'cohesion')
mfi = z_prop(pnt,'friction')*degrad
jc = z_prop(pnt,'jcohesion')
jfi = z_prop(pnt,'jfriction')*degrad
j_ang = beta * degrad ;;;;;;;;;degrad 常规变量 =pai/180
sm = 2.0 * mc * cos(mfi) / (1.0 - sin(mfi))
if betaj = 0 then
sj = -1.
else
sj = 2.0*jc / ((1.0-tan(jfi)*tan(j_ang))*sin(2.0*j_ang))
end_if
if sj < 0 then
anasol = sm
else
anasol = min(sj,sm)
end_if
end
;
def sigmav
sum=0.
pnt = gp_head
n = 0
loop while pnt # null
if gp_ypos(pnt) < eps then
sum = sum + gp_yfunbal(pnt)
n = n + 1
end_if
pnt = gp_next(pnt)
end_loop
sum = sum / c_area
sigmav = -sum
p_sol = anasol
p_err = -100. * (sum + p_sol) / p_sol
end
;
def servo
while_stepping
if step > 4500 * kk + 200
if unbal > high_unbal then
pnt = gp_head
loop while pnt # null
if gp_ypos(pnt) < eps then
gp_yvel(pnt) = gp_yvel(pnt) * 0.975
end_if
if gp_ypos(pnt) > 15.0 - eps then
gp_yvel(pnt) = gp_yvel(pnt) * 0.975
end_if
pnt = gp_next(pnt)
end_loop
end_if
if unbal < low_unbal then
pnt = gp_head
loop while pnt # null
if gp_ypos(pnt) < eps then
gp_yvel(pnt) = gp_yvel(pnt) * 1.025
if abs(gp_yvel(pnt)) > high_vel then
gp_yvel(pnt) = high_vel
end_if
end_if
if gp_ypos(pnt) > 15.0 - eps then
gp_yvel(pnt) = gp_yvel(pnt) * 1.025
if abs(gp_yvel(pnt)) > high_vel then
gp_yvel(pnt) = - high_vel
end_if
end_if
pnt = gp_next(pnt)
end_loop
end_if
end_if
end
;
set high_unbal 4 low_unbal 1
set high_vel 1.e-7
hist n 4500
hist sigmav
hist p_sol
hist beta
hist gp zdisp 0 0 0
hist p_err
test0
trisol
save triax.sav
ret
|
|