任意图形几何性质的数值算法

　　摘要：求解截面几何性质是一般结构分析的基础，对于普通的简单截面可以用解析方法求解，而对于任意复杂截面则必须采用数值方法，本文探讨了一种基于三角形单元的截面性质求解方法。
　　关键词：截面几何性质，数值方法
　　1引言
　　从材料力学的知识可知，截面几何性质的解析算法为积分形式，除了极少数规则形状截面，如圆形、矩形可简化成多项式形式以外，大部分非规则形状截面的几何性质都难以甚至不可能简化成多项式形式，这给应用带来了不小的麻烦，所以寻找数值方法来计算截面几何性质变得非常现实和必要。数值方法的基本原理就是先将任意形状的几何图形用多边形拟合，多边形的边数根据精度需要确定，再将此多边形离散成小的规则几何图形如是矩形单元（如图1.1）或者正方形单元（如图1.2），然后对每个小单元分别计算截面几何性质，最后通过移轴定理，将每个小单元的截面性质叠加起来，最终得出整个截面的几何性质。但是把任意图形离散成三角形单元则能以较小的单元数量得到较高的计算精度，本文推导了将任意截面离散成三角形的几何性质算法公式和数值算法。
　　2矩形单元和正方形单元的缺点
　　图2.1是采用矩形单元来计算截面几何性质的原理图，即将截面分割成宽度相等（如图中的S）的长条，然后将每个长条近似成矩形，这种方法要得到比较满意的精确结果，必须将单元分得足够细（即矩形宽度要小），否则会带来不小的误差，比如图2.1中的A、C两个单元，形状和标准的矩形相差甚远，如果要强制把它们近似成矩形，就会带来不小的误差。另外A、C单元刚好原离形心，所以它们对整个截面的几何性质影响也比在中间的B单元要大，但是却得不到比较精确的解。
　　
　　图2.1矩形单元法图2.2正方形单元法
　　图2.2是采用正方形单元来计算截面几何性质的原理图，显然，在细分成同样宽度的情况下，该方法精度是高于矩形单元法的，但是，代价就是平方倍的增加了计算量。对于图2.1计算量仅27个单元，而对于图2.2计算量却达到521个，计算量大大增加。另外，对于位于几何图形边缘的单元，如D、E单元近似为正方形也会带来误差，当然，可以采取一定措施来消除这种误差，比如超过1/2正方形的单元就当成正方形（如D单元），不超过1/2正方形的单元就忽略不计（如E单元）。这样正负相抵，误差大大减小。除非一些极端情况，位于边缘的单元全部是超过1/2正方形的，或者是全部都没有超过1/2正方形的，所以这种方法有时候也会带来不小的误差。
　　为了提高数值计算的精度，不得不更进一步的细分单元，对于复杂图形（如图2.1或2.2中的那种图形），这种做法倒不会太浪费，但如果对于简单图形（如图2.3或2.4那种图形），过多细分单元显得有点浪费。如图2.3这样的五边形，如果分成29个矩形单元，可得面积为36819.76，如果按右边分成413个单元，可得面积为36900，而该多边形的精确面积为36980.48，两种方法的误差分别为-0.43%和-0.22%。相对于采用矩形单元，采用正方形单元，单元数量增加了14.2倍，精度只提高了0.21%。详其它参数详见表2.1。
　　表2.1
　　 采用矩形单元时的参数 采用正方形单元的参数 精确的参数 采用矩形单元的误差 采用正方形单元的误差
　　单元数量 29 413 -- -- --
　　面积(m2) 36819.76 36900 36980.48 -0.43% -0.22%
　　Wx(m3) 730535.54 715152 712979 2.46% 0.30%
　　Wy(m3) 1020968.58 1022214.16 1026130.54 -0.50% -0.38%
　　Ix(m4) 88164573.65 87552920.06 88303989.31 -0.16% -0.85%
　　Iy(m4) 143595970.93 147168367.25 147049939.81 -2.35% 0.08%
　　
　　
　　图2.3矩形单元法（29个单元）图2.4正方形单元法（413个单元）
　　那有没有更好的单元划分方法呢？下面就讨论一种基于三角形单元划分的数值方法。
　　3三角形单元算法的原理
　　正如上述讨论的，采用矩形单元或者是正方形单元，要么单元过多，计算量偏大，要么精度不高，其原因就是因为我们采用的单元只是多边形的近似。如果将多边形离散成三角形单元，则可以大大减少这种因近似带来的误差。如图3.1所示，将多边形离散成三角形单元，不仅单元数量大大减少，而且计算精度也大为提高，根据图3.1划分的单元计算得到的多边形面积为36980.48，误差为0.0%。
　　
　　图3.1三角形单元（5个单元）图3.2形心
　　
　　3.1材料力学公式
　　常用的截面几何性质有如下几种：
　　3.1.1形心（如图3.2）
　　公式如下：
　　
　　3.1.2静矩
　　公式如下：
　　
　　3.1.3惯性矩
　　公式如下：
　　
　　
　　
　　3.2三角形单元截面性质的推导
　　推导之前，做如下假定，三角形的一个顶点位于坐标原点（如图3.4），至于不在坐标原点的情况（如图3.3），可以通过坐标变换，将其转换到坐标原点上。假设XOY平面上任意三角形三个顶点的坐标分别为：O'（xo,yo)，A(xa,ya)，B(xb,yb)。我们规定，xo≤xa≤xb，现通过坐标变换，将O'点作为新坐标系的原点，则得到此三角形在新坐标系xO'y下的坐标O'(0,0)，A（xa-xo，ya-yo），B（xb-xo,yb-yo）。
　　
　　图3.3整体坐标系下的三角形图3.4局部坐标系下的三角形
　　经推导，得到如下关于三角形局部坐标轴xO'y的截面参数：
　　（1）面积：         ……（公式1）
　　（2）对于x轴的静矩:      ……（公式2）
　　                 ……（公式3）
　　              ……（公式4）
　　（3）对于y轴的静矩:       ……（公式5）
　　  ……（公式6）
　　              ……（公式7）
　　（4）形心：             ……（公式8）
　　               ……（公式9）
　　（5）对于x轴惯性矩：      ……（公式10）
　　
　　 ……（公式11）
　　              ……（公式12）
　　（6）对于y轴惯性矩：       ……（公式13）
　　   ……（公式14）
　　              ……（公式15）
　　以上公式中xa，xb，xc，ya，yb，yc均为坐标，ka，kb，kab均为斜率，其中：
　　ka=ya/xa
　　kb=yb/xb
　　kab=(ya-yb)/(xa-xb)
　　3.3任意图形截面性质的推导
　　对于任意复杂图形的截面性质，只需要将各个小三角形单元的几何性质“叠加”起来。但是不是普通的相加，方法如下：
　　（1）面积：……（公式16），其中n为三角形单元总数，Ai为各个三角形单元面积，按公式1计算；
　　（2）形心:         ……（公式17）
　　……（公式18），其中n为三角形单元总数，Ai为各个三角形单元面积，按公式1计算；
　　（3）对于x轴的静矩:         ……（公式19）
　　（4）对于y轴的静矩:         ……（公式20）
　　（5）对于x轴惯性矩：……（公式21），其中n为三角形单元总数，Ixi为各个三角形单元对于x轴的惯性矩，按公式10~12计算；
　　（6） 对于y轴惯性矩：……（公式22），其中n为三角形单元总数，Iyi为各个三角形单元对于y轴的惯性矩，按公式13~15计算；
　　4结语
　　本文推导了任意几何图形划分成三角形单元时的截面性质公式。并比较了这种方法相对于其他方法的优缺点。本方法算法简单、快速、精确，容易编制成计算机程序，方便使用。
　　参考文献：
　　[1]苏翼林主编《材料力学》天津大学出版社2001.6
　　[2]同济大学数学教研室主编《高等数学》高等教育出版社第四版
　　[3]郝际平钟炜辉《薄壁杆件的弯曲与扭转》高等教育出版社2006.9