[讨论]：给各位斑竹出难题;能不能编个计算极惯性矩的LISP

MJJ3468

这个太难了，好像目前都没有（有没做结构的啊）
非圆截面的极惯性矩（抗扭惯性矩）
目前也只能近似的计算，而且非常的麻烦，除非用有限元软件。
那位斑竹能这方面的程序？？哈哈。这个还真让人闹心

xyp1964

不知可行：

2. [FONT=courier new];;;封闭域的形心、惯性矩、惯性积、旋转半径
   
3. (defun c:test285 ()
   
4.   (defun write-tmp ()
   
5.     (setq f (open "c:/tmp.mpr" "r"))
   
6.     (read-line f)
   
7.     (read-line f)
   
8.     (read-line f)
   
9.     (while (setq txt (read-line f))
   
10.       (setq pt (list (car pt) (- (cadr pt) (* sc 800))))
    
11.       (xyp-Text 1 pt txt)
    
12.     )
    
13.     (close f)
    
14.   )
    
15.   (cmdlasc0)
    
16.   (setvar "osmode" 0)
    
17.   (xyp-mkLaCo "中心线" 1)
    
18.   (while (setq pt (getpoint "\n选择区域内一点<退出> : "))
    
19.     (command "-boundary" "a" "o" "r" "" pt "")
    
20.     (setq s1 (entlast)
    
21.           pt (vlax-safearray->list
    
22.                (vlax-variant-value
    
23.                  (vla-get-centroid (xyp-Ve2Vo s1))
    
24.                )
    
25.              )
    
26.     )
    
27.     (xyp-Cross pt (* sc 200) 1)
    
28.     (command "MASSPROP" s1 "" "y" "c:/tmp.mpr")
    
29.     (xyp-mkLaCo "面域标注" 4)
    
30.     (write-tmp)
    
31.     (redraw)
    
32.   )
    
33.   (cmdla1)
    
34. )[/FONT]

[FONT=courier new];;;封闭域的形心、惯性矩、惯性积、旋转半径<br /> (defun c:test285 ()<br />   (defun write-tmp ()<br />     (setq f (open "c:/tmp.mpr" "r"))<br />     (read-line f)<br />     (read-line f)<br />     (read-line f)<br />     (while (setq txt (read-line f))<br />       (setq pt (list (car pt) (- (cadr pt) (* sc 800))))<br />       (xyp-Text 1 pt txt)<br />     )<br />     (close f)<br />   )<br />   (cmdlasc0)<br />   (setvar "osmode" 0)<br />   (xyp-mkLaCo "中心线" 1)<br />   (while (setq pt (getpoint "\n选择区域内一点<退出> : "))<br />     (command "-boundary" "a" "o" "r" "" pt "")<br />     (setq s1 (entlast)<br />           pt (vlax-safearray->list<br />                (vlax-variant-value<br />                  (vla-get-centroid (xyp-Ve2Vo s1))<br />                )<br />              )<br />     )<br />     (xyp-Cross pt (* sc 200) 1)<br />     (command "MASSPROP" s1 "" "y" "c:/tmp.mpr")<br />     (xyp-mkLaCo "面域标注" 4)<br />     (write-tmp)<br />     (redraw)<br />   )<br />   (cmdla1)<br /> )[/FONT]

MJJ3468

最初由 xyp1964 发布
[B]不知可行：
[CODE];;;封闭域的形心、惯性矩、惯性积、旋转半径
(defun c:test285 ()
  (defun write-tmp ()
    (setq f (open "c:/tmp.mpr" "r"))
    (read-line f)
    (read-line f)
... [/B]

很遗憾，不行，你这就是CAD自带的查属性的。但不能计算极惯性矩也就是常说的IP。
众所周知，EA为拉压刚度，EI为抗弯刚度，GJ为抗扭刚度，SaG为抗剪刚度，除了抗扭J,其他都可以通过CAD啊，和其他的截面特性程序方便查询，J其实就是截面的极惯性矩，圆的J=IP=Ix+Iy,其他典型的截面有相应的公式，但对于不规则的就无从下手，请问论坛大虾有啥办法，这个是太难了

飞马牌汽车

看看这个吧
http://www.slickwin.com/massplus.html
是lisp的，虽是共享，可以用
把http://www.slickwin.com/mpftp/mpwin.zip
解压到和acad.exe同目录下，先用mp，再用answer
其中
Jz         Polar moment of Inertia about the centroid

MJJ3468

最初由 飞马牌汽车 发布
[B]看看这个吧
http://www.slickwin.com/massplus.html
是lisp的，虽是共享，可以用
把http://www.slickwin.com/mpftp/mpwin.zip
解压到和acad.exe同目录下，先用mp，再用answer
其中
Jz  ... [/B]

楼上的大虾，这个都被你找到，真牛。
不过，要许可，才能计算正确，请看下面介绍（我用了下，计算的全不对，连面积都不正确） ：P

飞马牌汽车

这个软件的限制似乎不在这里，而在于剖分的精度，起码我这边是这样的
极惯性矩确实要进行剖分计算的，lisp应该还是要编好一段的
既然是结构方面的，不如尝试ansys, sap2000等大软件，其内置的功能都能计算
或者csi的section builder或者其他的shapebuilder
section builder做这个可以，如图，软件请自己找找

MJJ3468

最初由 飞马牌汽车 发布
[B]这个软件的限制似乎不在这里，而在于剖分的精度，起码我这边是这样的
极惯性矩确实要进行剖分计算的，lisp应该还是要编好一段的
既然是结构方面的，不如尝试ansys, sap2000等大软件，其内置的功能都能计算
或者cs... [/B]

ansys,sap2000,各版本的我都有，我知道都能算，但不是为了更简便一点吗，为了查惯性矩，谁会去用ANSYS呢？哈哈，谢谢。

飞马牌汽车

:) 各版本都有就太吓人了
利用以前的一些Lisp程序改造了一下
有一点积分技术，采用1000根左右扫描线，速度颇慢，要7，8秒吧，针对region封闭无孤岛物体进行。由于程序结构乱七八糟的，就不献丑了，提供vlx测试一下吧，反正将就能用。
由于不懂采用双精度计算，因此总有一点点误差，有空再修改了。

MJJ3468

最初由 飞马牌汽车 发布
[B]:) 各版本都有就太吓人了
利用以前的一些Lisp程序改造了一下
有一点积分技术，采用1000根左右扫描线，速度颇慢，要7，8秒吧，针对region封闭无孤岛物体进行。由于程序结构乱七八糟的，就不献丑了，提供vlx测试一下... [/B]

原来是华南理工的才子，不知你是否是学结构的。
测试了下，不是一点点误差，是很大。哈哈

飞马牌汽车

不知道您的误差如何，可以告知我么
如有数据，请告知我，我测试一下，看是否哪里有错误，自当改进，请别嘲笑。
Ix是对形心的，不是对x轴的，这点我忘记说了
我做的直径100的圆的测试如下
这个误差与双精度计算有关，与扫描精度取值也有关系。

MJJ3468

只有圆截面IP=IX+IY
下面已H型钢为例。
模型为175X175X7.5X11
你的程序计算结果
ANSYS的计算结果，TORSION CONSTANT为扭转常量，即极惯性矩

snoopychen

谢谢你：）
我重新思考了极惯性矩这个问题，在lisp版还是不想用“汽车”的马甲了，还是会食言
我一直以来认为Ip＝积分（ρ^2dA），按道理来说一般比两个方向的惯性矩都应该大些
材料力学书上也是这样定义的，孙训方《材料力学》P104。
因而程序计算的时候是按照如下步骤进行：
1.找出形心
2.用水平线扫描region，得到2个或者更多个交点，任意两个交点之间求出一条线段，对于求面积

的时候，取出这个线段的长度，乘以每次的增量高度，得到微元面积。对于惯性矩，通过对面积元

与形心y方向距离相乘，叠加而到。而对于极惯性矩，采用了对线段再次分为微元，依次寻找微元与

形心的距离积分，得到这个结果。应该说结果是按照Ip＝积分（ρ^2dA）来得到的。

但是就如你所说的，扭转常数却经常和这个数值不同。除了圆形之外，可能得用其他的方法来计算

这个了。
比如下图1的方形，其极惯性矩可以用斜段ρ来计算，但是在进行抗扭计算的时候，其剪应力的方向

其实不是垂直ρ的（那个蓝色箭头），而是沿着水平方向（黄色箭头），而对于薄壁杆件，如图3的

工字形截面，其剪应力更是沿着曲线的切线方向走。所以在你刚才提出ansys的结果不同之后，我

用section builder重新计算之后，发现其工字形截面扭转常量小于普通惯性矩之后，开始还有点想

不明白。其实它的计算涉及到了剪应力的方向。
非常谢谢你，我现在觉得极惯性矩和扭转常量之间可能不是简单的关系了，这个问题我一时也解决

不了，等我查查书，想明白了再看看能否编写。
有空多探讨。下面是初始代码，请指正。
[php]
;;some of the following code are written by qjchen
;;purpose: to calculate the Ip, polar moment inertia
(defun c:test (/ disall iall ijall ent cena p1 cenx ceny a b dis ay by th th0 lp
rp inter1 i len j dis12 area disy iall dis34 dis56 ijall)
(vl-load-com)
(setq ijall 0 iall 0 area 0)
  (startTimer)
  (setting)
  (prompt "\n请选择面域region:" )
  (setq ent (car (entsel)))
  (setq cena (cent ent))
  (setq p1 (xyval1 ent))
  (setq cenx (nth 0 cena))
  (setq ceny (nth 1 cena))  
  (setq a (polar (nth 0 p1) (* 0.75 pi) 1))
  (setq b (polar (nth 1 p1) (* 1.75 pi) 1))  
  (setq dis (/ (distance a b) 1000))
  (setq ay (nth 1 a)
        by (nth 1 b)
  )
  (setq th by)
  (setq th0 dis)
  (while (< th ay)
    (setq lp (list (nth 0 a) th 0))
    (setq rp (list (nth 0 b) th 0))
    ;(grdraw lp rp 249)
    (setq inter1 (vl-Get-Int-Pt lp rp))
   (setq i 0
          len (length inter1)
    )
    (if (> len 0)
      (if (/= (* (/ len 2) 2) len)
        (setq inter1 (append
                       inter1
                       (list (last inter1))
                     )
        )
      )
    )
    (setq len (length inter1))
    (setq j 0)
    (repeat (/ len 2)
      (setq dis12 (* (distance (nth j inter1) (nth (1+ j) inter1)) dis))
      (setq area (+ area dis12))
      (setq disy (abs (- th ceny)))
      (setq iall (+ iall (* dis12 (* disy disy))))
      (setq dis34 (distance (nth j inter1) (nth (1+ j) inter1)))
      (setq dis56 (abs (- (nth 0 (nth j inter1)) cenx)))
      (setq ijall (+ ijall (int2 dis34 dis56 disy dis)))
      (setq j (+ j 2))
    )
    (setq th (+ th th0))
  )
  (resetting)
  (endTimer (vl-symbol-name 'c:test))
(princ "\n 面积是：") (princ area);Erase_DV
(princ "\n 绕x轴形心的惯性矩：") (princ iall);Erase_DV
(princ "\n 绕形心的极惯性矩：") (princ ijall);Erase_DV
(princ)
)

;;;calculate a line moment inertia to a point
(defun int2 (a b c tt / f1 f2 f3 f4 f5)
  (setq f1 (* c c a))
  (setq f2 (* b b a))
  (setq f3 (* b a a))
  (setq f4 (/ (* a a a) 3))
  (setq f5 (- (+ f1 f2 f4) f3))
  (setq f5 (* f5 tt))
  f5
)
;;cal the centroid of the selected region
(defun cent (ent / b c centroid i d)
  (setq b (vlax-ename->vla-object ent))
  (setq centroid (vla-get-centroid b))
  (setq i 0
        d nil
  )
  (repeat 2
    (setq d (append
              d
              (list (vlax-safearray-get-element (vlax-variant-value centroid)
                                                i
                    )
              )
            )
    )
    (setq i (1+ i))
  )
  d
)
;;; -------------------------------------------------------------------
;;; | The following code are taken from xoutside.com                  |
;;; | http://xoutside.com/CAD/lisp/lisp_chair.htm                     |
;;; | Thanks to the Korea friend                                      |
;;; | Purpose: Get the intersection of Two object                     |
;;; -------------------------------------------------------------------
(defun vl-Get-Int-Pt (FirstPoint SecondPoint / acadDocument mSpace SSetName
                                 SSets SSet reapp ex obj Baseline
                     )
  (vl-load-com)
  (setq acadDocument (vla-get-ActiveDocument (vlax-get-acad-object)))
  (setq mSpace (vla-get-ModelSpace acadDocument))
  (setq SSetName "MySSet")
  (setq SSets (vla-get-SelectionSets acadDocument))
  (if (vl-catch-all-error-p (vl-catch-all-apply 'vla-add (list SSets
                                                               SSetName
                                                         )
                            )
      )
    (vla-clear (vla-Item SSets SSetName))
  )
  (setq SSet (vla-Item SSets SSetName))
  (setq Baseline (vla-Addline mspace (vlax-3d-point FirstPoint)
                              (vlax-3d-point SecondPoint)
                 )
  )
  (vla-SelectByPolygon SSet acSelectionSetFence
                       (kht:list->safearray (append
                                              FirstPoint
                                              SecondPoint
                                            ) 'vlax-vbdouble
                       )
  )
  (vlax-for obj sset (if (setq ex (kht-intersect
                                                 (vlax-vla-object->ename BaseLine)
                                                 (vlax-vla-object->ename obj)
                                  )
                         )
                       (setq reapp (append
                                     reapp
                                     ex
                                   )
                       )
                     )
  )
  (vla-delete BaseLine)
  (setq reapp (vl-sort reapp '(lambda (e1 e2)
                                (< (car e1) (car e2))
                              )
              )
  )
  reapp
)
;;; Change to avoid text intersection
(defun kht-intersect (en1 en2 / a b x ex ex-app c d e)
  (vl-load-com)
  (setq c (cdr (assoc 0 (entget en1)))
        d (cdr (assoc 0 (entget en2)))
  )
  (if (or
        (= c "TEXT")
        (= d "TEXT")
      )
    (setq e -1)
  )
  (setq En1 (vlax-ename->vla-object En1))
  (setq En2 (vlax-ename->vla-object En2))
  (setq a (vla-intersectwith en1 en2 acExtendNone))
  (setq a (vlax-variant-value a))
  (setq b (vlax-safearray-get-u-bound a 1))
  (if (= e -1)
    (setq b e)
  )
  (if (/= b -1)
    (progn
      (setq a (vlax-safearray->list a))
      (repeat (/ (length a) 3)
        (setq ex-app (append
                       ex-app
                       (list (list (car a) (cadr a) (caddr a)))
                     )
        )
        (setq a (cdr (cdr (cdr a))))
      )
      ex-app
    )
    nil
  )
)
(defun kht:list->safearray (lst datatype)
  (vlax-safearray-fill (vlax-make-safearray (eval datatype) (cons 0
                                                                  (1-
                                                                      (length lst)
                                                                  )
                                                            )
                       ) lst
  )
)
;;; ----------------------------------------------------------
;;; |           midpoint function                            |
;;; ----------------------------------------------------------
(defun midp (p1 p2)
  (mapcar
    '(lambda (x)
       (/ x 2.)
     )
    (mapcar
      '+
      p1
      p2
    )
  )
)
;;; -----------------------------------------------------------------
;;; | The following code taken from Mr.Tony Hotchkiss at Cadalyst   |
;;; | To set and reset the system variable                          |
;;; -----------------------------------------------------------------
(defun err (s)
  (if (= s "Function cancelled")
    (princ "\nfunction - cancelled: ")
    (progn
      (princ "\nfunction - Error: ")
      (princ s)
      (terpri)
    )        ; _ end of progn
  )        ; _ end of if
  (resetting)
  (princ "SYSTEM VARIABLES have been reset\n")
  (princ)
)
;;; err
;;; setting and resetting the system variables
(defun setv (systvar newval / x)
  (setq x (read (strcat systvar "1")))
  (set x (getvar systvar))
  (setvar systvar newval)
)
;;; setv
(defun setting ()
(command "_undo" "_be")
  (setq oerr *error*)
  (setq *error* err)
  (setv "BLIPMODE" 0)
  (setv "CMDECHO" 0)
  (setv "OSMODE" 0)
)
;;; setting
(defun rsetv (systvar)
  (setq x (read (strcat systvar "1")))
  (setvar systvar (eval x))
)
;;; rsetv
(defun resetting ()
  (rsetv "BLIPMODE")
  (rsetv "CMDECHO")
  (rsetv "OSMODE")
  (setq *error* oerr)
  (command "_undo" "_e")
)
;;; -----------------------------------------------------------------
;;; | The following code taken from www.theswamp.org                |
;;; | To calculate the time that the program run                    |
;;; -----------------------------------------------------------------
(defun startTimer ()
  (setq time (getvar "DATE"))
)
(defun endTimer (func)
  (setq time (- (getvar "DATE") time)
        seconds (* 86400.0 (- time (fix time)))
  )
  (gc)
  (outPut seconds func)
)
(defun outPut (secs def)
  (princ "\nProgressing...")
  (gc)
  (princ (strcat "\n 总耗时: " (rtos secs 2 6) "秒"))
  (princ)
)
;(princ "\n Please use the test command to run")
;;; _ end of xyval                                                   ;
;;;;---The following codes are copy from Tony Hotchkiss at cadalyst  ;
;;Get the boundingbox of one object                                  ;
(defun xyval1 (obj / minpt maxpt topy bottmy leftx rightx)
  (setq obj (vlax-ename->vla-object obj))
  (vla-GetBoundingBox obj 'minpt 'maxpt)
  (setq pt1 (vlax-safearray->list minpt)
        pt2 (vlax-safearray->list maxpt)
        topy (cadr pt2)
        bottmy (cadr pt1)
        leftx (car pt1)
        rightx (car pt2)
  )                                       
  (list (list leftx topy) (list rightx bottmy))
)


(princ "\n查询极惯性矩，命令test")
(princ "\n只能对region面域进行查询")
(princ)

[/php]

and2008

很需要这方面的程序哈，期待早日完善..