求助英文手册里的单轴压缩的例子看不懂，求高手详解一下

lzh423343315

就是英文手册后面的一个例子。主要是里面定义的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

lzh423343315

求高手给解释一下，详细点，我这初学，自己琢磨，有些东西确实书上也找不到，手册里也查不明白。谢谢啊！

lzh423343315

给点解释，谢谢了。。。

lzh423343315

我的个乖乖。。。。没人会么？

lzh423343315

一天了还没人答复

legealdragon

前三个是定义的函数
后面的是设置最大不平衡力、最小不平衡力和最大节点速度

lzh423343315

谢谢了。还在继续研究。

weason

貌似十分复杂，这个例子，楼主看懂了吗？