本帖最后由 745257116 于 2015-12-29 10:52 编辑
声明:此命令流为米兰的小铁匠原创,也谢谢作者在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 |