FLAC实体单元求内力（弯矩、轴力、剪力）

bob410

FLAC实体单元（ZONE）,求截面内力（弯矩、轴力、剪力），有好的方法没？

289869308

同求啊……

东701

书上应该有

machi

声明：此命令流为米兰的小铁匠原创，也谢谢作者在QQ空间上的分享。
以下为原版内容：
; 存档
; Flac3D合成环向弯矩和轴力


demo.txt

call split_liner.fis
range name primary group 2
range name outer group 1
range name inner group 3
set @grname = 'primary'
set @outergr = 'outer'
set @innergr = 'inner'
@split_liner
call get_mn_annular.fis
set @mn_file = 'flac3d_mn_primary.csv'
set @sec_th = 0.28
set @elem_le= 1.0
set @xx_cen = 0.0
set @zz_cen = 0.0
set @ksi = 0.5
range name flac3dmn_range1 group primary_outer
range name flac3dmn_range2 group primary_inner
@get_mn_annular




split_liner.fis

def split_liner
;grname = 'primary'
;outergr = 'outer'
;innergr = 'inner'
p_z = zone_head
loop while p_z # null
flag1 = 0
flag2 = 0
if in_range(grname, p_z) = 1 then
gpcount = z_numgp(p_z)
zid = z_id(p_z)
loop i(1, gpcount)
p_gp = z_gp(p_z, i)
flag1 = flag1 + in_range(outergr, p_gp)
flag2 = flag2 + in_range(innergr, p_gp)
endloop
section
flag3 = flag1 * flag2
if flag3 # 0 then
exit section
endif
if flag1 > 0 then
newgrname = grname + '_' + outergr
z_group(p_z) = newgrname
exit section
endif
if flag2 > 0 then
newgrname = grname + '_' + innergr
z_group(p_z) = newgrname
endif
endsection
endif
p_z = z_next(p_z)
endloop
end


get_mn_annular.fis

def get_mn_annular
; number of zones, address
n_zone1 = 0
n_zone2 = 0
m_pnt1 = null
m_pnt2 = null
z_pnts = z_list
loop foreach z_pnt z_pnts
if in_range('flac3dmn_range1', z_pnt) = 1 then
n_zone1 = n_zone1 + 1
m_pnt11 = get_mem(2)
mem(m_pnt11) = m_pnt1
mem(m_pnt11+1) = z_pnt
m_pnt1 = m_pnt11
endif
if in_range('flac3dmn_range2', z_pnt) = 1 then
n_zone2 = n_zone2 + 1
m_pnt21 = get_mem(2)
mem(m_pnt21) = m_pnt2
mem(m_pnt21+1) = z_pnt
m_pnt2 = m_pnt21
endif
endloop
; open
status = open(mn_file, 1, 1)
buf = get_array(1)
buf(1) = string(xx_cen)+','+string(zz_cen)
status = write(buf, 1)
buf(1) = 'x0,z0,x1,z1,x2,z2,m,n'
status = write(buf, 1)
; resultant moment and nforce
loop while m_pnt11 # null
p_z = mem(m_pnt11 + 1)
x_cor1 = z_xcen(p_z)
z_cor1 = z_zcen(p_z)
sgm_xx1= z_sxx(p_z)
sgm_xz1= z_sxz(p_z)
sgm_zz1= z_szz(p_z)
; mpnt21
m_pnt2 = m_pnt21
p_zz = mem(m_pnt2 + 1)
x_cor2 = z_xcen(p_zz)
z_cor2 = z_zcen(p_zz)
sgm_xx2= z_sxx(p_zz)
sgm_xz2= z_sxz(p_zz)
sgm_zz2= z_szz(p_zz)
sd_dist= sqrt((x_cor1-x_cor2)^2 + (z_cor1-z_cor2)^2)
p_z2 = p_zz
loop while m_pnt2 # null
p_zz = mem(m_pnt2 + 1)
x_cor3 = z_xcen(p_zz)
z_cor3 = z_zcen(p_zz)
sd_dist3= sqrt((x_cor1-x_cor3)^2 + (z_cor1-z_cor3)^2)
if sd_dist > sd_dist3 then
p_z2 = p_zz
sd_dist = sd_dist3
x_cor2 = x_cor3
z_cor2 = z_cor3
sgm_xx2 = z_sxx(p_zz)
sgm_xz2 = z_sxz(p_zz)
sgm_zz2 = z_szz(p_zz)
endif
m_pnt2 = mem(m_pnt2)
endloop
; get elem_length
gps = get_array(8, 3)
loop ii(1, 8)
gps(ii, 1) = z_gp(p_z, ii)
gps(ii, 2) = z_gp(p_z2,ii)
gps(ii, 3) = 0
endloop
loop ii(1, 8)
loop iii(1,8)
if gps(ii,1) = gps(iii,2) then
gps(ii, 3) = 1
endif
endloop
endloop
iii = 1
_gps = get_array(2)
loop ii(1, 8)
if gps(ii, 3) = 1 then
if gp_ypos(gps(ii,1)) < z_ycen(p_z) then
gps(ii, 3) = 2
_gps(iii) = gps(ii, 1)
iii = iii + 1
endif
endif
endloop
delta_x = gp_xpos(_gps(1))-gp_xpos(_gps(2))
delta_z = gp_zpos(_gps(1))-gp_zpos(_gps(2))
; elem_le = sqrt(delta_x ^ 2 + delta_z ^ 2)
ii = lose_array(gps)
ii = lose_array(_gps)
sec_th = 2 * sd_dist
; θ = arctan(dx/dy)
alpha=-atan((x_cor1-x_cor2)/(z_cor1-z_cor2))
; σn = σx*cos(θ)^2 + σy*cos(θ)^2 + σxy*sin(2θ)
sgm_n1=sgm_xx1*cos(alpha)^2 + sgm_zz1*sin(alpha)^2 + sgm_xz1*sin(2*alpha)
sgm_n2=sgm_xx2*cos(alpha)^2 + sgm_zz2*sin(alpha)^2 + sgm_xz2*sin(2*alpha)
; σ1 = (σn1 + σn2) / 2 + (σn1 + σn2) / 2τ
sgm_sg1=0.5*(sgm_n1+sgm_n2)+(sgm_n1-sgm_n2)/(2.0*ksi)
; σ2 = (σn1 + σn2) / 2 - (σn1 + σn2) / 2τ
sgm_sg2=0.5*(sgm_n1+sgm_n2)-(sgm_n1-sgm_n2)/(2.0*ksi)
; m = (σ1 - σ2) * bh^2 / 12
m_moment=elem_le*sec_th*sec_th*(sgm_sg1-sgm_sg2)/12
; n = (σ1 + σ2) * bh / 2
n_force =elem_le*sec_th*(sgm_sg1+sgm_sg2)/2
sec_xcor=(x_cor1+x_cor2)/2
sec_zcor=(z_cor1+z_cor2)/2
buf(1)=buf(1)+string(sec_xcor)+','
buf(1)=buf(1)+string(sec_zcor)+','
buf(1)=buf(1)+string(x_cor1)+','
buf(1)=buf(1)+string(z_cor1)+','
buf(1)=buf(1)+string(x_cor2)+','
buf(1)=buf(1)+string(z_cor2)+','
buf(1)=buf(1)+string(m_moment)+','
buf(1)=buf(1)+string(n_force)+','
status = write(buf, 1)
; release memory
m_pnt1 = m_pnt11
m_pnt11 = mem(m_pnt11)
ii = lose_mem(2, m_pnt1)
endloop
; release memory
loop while m_pnt21 # null
m_pnt2 = m_pnt21
m_pnt21 = mem(m_pnt21)
ii = lose_mem(2, m_pnt2)
endloop
; close
status = close
ii = out('Successfully Write ' + string(n_zone1) + ' Data Into File ' + mn_file)
ii = lose_array(buf)
end
; @get_mn_annular
原帖地址
http://forum.simwe.com/forum.php ... =%E5%86%85%E5%8A%9B

1771864883

111rrrrrrrrrrrrrrrrrrrrrrrrr