面向ChatGPT编程：有限单元法求解任意截面矩
“ 由于有限单元法可以实现将任意截面划分为三角形单元，所以可以实现对任意截面的极惯性矩计算。本文将使用HyperMesh进行网格划分，并使用tcl语言进行二次开发，计算任意截面极惯性矩。由于本人对tcl语法不熟悉，对于不熟悉的知识将请教ChatGPT，以加快开发速度。”
为验证算法的正确性，先从简单的方钢开始，使用HyperMesh建立一个外宽60mm，外高40mm，壁厚5mm，使用三角形（1mm）网格完成网格划分。如下图所示：
解决问题的思路：
使用HyperMesh获得目标截面的离散化网格，使用三角形单元划分，因为三角形单元面积以及形心非常好计算，所以我们需要一些工具函数：
需要一个用于根据三个节点的坐标计算三角形面积的tcl函数；
需要一个函数根据两点坐标计算两点距离，用于确定单元距离形心的距离
我们的目标是计算惯性矩和极惯性矩，都需要用到整体截面形心的信息，所以需要先通过截面静矩计算得到截面形心坐标，然后根据微积分公式计算得到想要的目标。
01
—
tcl工具函数的准备
由于对tcl语言不熟悉的缘故，想要自行编写一些tcl工具函数是非常耗时的，基本上所有的命令都需要从手册或者网上搜索，，所以我们请教一下ChatGPT：
通过只能生成的程序需要检验其准确性，但是运行时HyperMesh提示没有vecsub函数，只能再问一下ChatGPT：
发现可能是版本的问题，并且他给出了函数实现方式，使用同样的方法获得其余的一些HyperMesh中不存在的函数，包括vecnorm、vecmag、vecprod：
获得了所有函数后，在此运行程序就没有报错了。现在需要确定其计算结果的正确性，经过检查其存在错误，计算过程不应该使用vecnorm进行归一化，修改此处错误后即可获得正确的面积以及两点距离，需要的工具函数如下，这些函数均是ChatGPT生成：
# 向量差proc vecsub {a b} {set result {}for {set i 0} {$i < [llength $a]} {incr i} {lappend result [expr {[lindex $b $i] - [lindex $a $i]}]}return $result}# 向量单位化proc vecnorm {a} {set mag [vecmag $a]if {$mag == 0} {return {0 0 0}}set result {}for {set i 0} {$i < [llength $a]} {incr i} {lappend result [expr {[lindex $a $i] / $mag}]}return $result}# 向量模长proc vecmag {a} {set sum 0foreach x $a {set sum [expr $sum + $x*$x]}return [expr sqrt($sum)]}# 向量叉积proc vecprod {a b} {set x [expr [lindex $a 1]*[lindex $b 2] - [lindex $a 2]*[lindex $b 1]]set y [expr [lindex $a 2]*[lindex $b 0] - [lindex $a 0]*[lindex $b 2]]set z [expr [lindex $a 0]*[lindex $b 1] - [lindex $a 1]*[lindex $b 0]]return [list $x $y $z]}
我们需要进行均值计算和求和计算，所以也要让ChatGPT帮我们写一下：
# 列表均值proc eva {mylist} {set sum 0foreach num $mylist {set sum [expr {$sum + $num}]}set mean [expr {$sum / [llength $mylist]}]return $mean}# 列表求和proc sum {a} {set result 0foreach element $a {set result [expr {$result + $element}]}return $result}
后续程序中需要多次根据节点ID获取节点坐标，所以这个功能我们也要包装一下，自行创建工具函数，函数结构可以借鉴ChatGPT生成的：
# 通过节点ID获取节点坐标proc get_node_coords {i} {set j [hm_nodevalue $i coordinates]set j [regsub -all "{" $j " "]set j [regsub -all "}" $j " "]return $j}
现在我们需要遍历选择出来的单元，这需要请教一下ChatGPT：
02
—
确定形心
静矩的定义：first moment of area。截面的静矩是指截面内每个点到某个轴线距离的平方与该点所对应的面积的乘积之和。公式为：
其中，A是截面的面积，y是每个点到x轴线的距离。
静矩作用：可用于计算截面形心。截面对某轴的静矩为零，则该轴必过形心，截面对一个坐标系的两个轴的静矩都为零，则该坐标系原点为形心。
# 选择截面单元*createmarkpanel elems 1 "Please select the elements of a surface."set es [hm_getmark elems 1]# 获取所有节点并去重set nodes {}# 获取所有节点坐标，并计算z坐标均值set zs {}# 对X轴静矩set SX {}# 对Y轴静矩set SY {}# 总面积set Area {}foreach e $es {set n1 [hm_getvalue elems id=$e dataname=node1]set A [get_node_coords $n1]set n2 [hm_getvalue elems id=$e dataname=node2]set B [get_node_coords $n2]set n3 [hm_getvalue elems id=$e dataname=node3]set C [get_node_coords $n3]lappend nodes $n1lappend nodes $n2lappend nodes $n3# puts "$A $B $C"# 计算三角形面积set AB [vecsub $B $A]set AC [vecsub $C $A]set cross [vecprod $AB $AC]set area [expr 0.5 * [vecmag $cross]]lappend Area $area# 计算单元形心坐标set Gx [expr ([lindex $A 0] + [lindex $B 0] + [lindex $C 0])/3.0]set Gy [expr ([lindex $A 1] + [lindex $B 1] + [lindex $C 1])/3.0]set Gz [expr ([lindex $A 2] + [lindex $B 2] + [lindex $C 2])/3.0]# 计算静矩，用于确定形心坐标set Sx [expr $Gy * $area]lappend SX $Sxset Sy [expr $Gx * $area]lappend SY $Sy}# 根据xy轴的静矩以及总面积计算形心xy坐标， 根据z轴坐标均值确定形心z轴坐标set nodes [lsort -unique $nodes]foreach n $nodes {set D [get_node_coords $n]lappend zs [lindex $D 2]}set staticmomentx [sum $SX]set staticmomenty [sum $SY]set Area [sum $Area]set x0 [expr $staticmomenty/$Area]set y0 [expr $staticmomentx/$Area]set z0 [eva $zs]puts "形心： $x0 $y0 $z0"
03
—
惯性矩和极惯性矩
惯性矩定义：area moment of inertia，截面上所有点至坐标轴距离平方的和，可反映截面上的点相对于轴的分布情况。
惯性矩作用：惯性矩可用于计算纯弯曲变形杆截面上的正应力
极惯性矩定义：polar moment of inertia，面上所有点到某点的距离平方的和，可反映截面上的点相对于某点的分布情况，计算公式为
极惯性矩作用：极惯性矩可用于计算杆件在扭转状态下的最大切应力。
# 极惯性矩set IZ {}# 对X轴惯性矩set IX {}# 对Y轴惯性矩set IY {}foreach e $es {set n1 [hm_getvalue elems id=$e dataname=node1]set A [get_node_coords $n1]set n2 [hm_getvalue elems id=$e dataname=node2]set B [get_node_coords $n2]set n3 [hm_getvalue elems id=$e dataname=node3]set C [get_node_coords $n3]# 计算三角形面积set AB [vecsub $B $A]set AC [vecsub $C $A]set cross [vecprod $AB $AC]set area [expr 0.5 * [vecmag $cross]]# puts $AB# puts $AC# puts $area# puts "顶点：$A $B $C"# 计算单元形心坐标set Gx [expr ([lindex $A 0] + [lindex $B 0] + [lindex $C 0])/3.0]set Gy [expr ([lindex $A 1] + [lindex $B 1] + [lindex $C 1])/3.0]set Gz [expr ([lindex $A 2] + [lindex $B 2] + [lindex $C 2])/3.0]set distance [vecmag [vecsub "$Gx $Gy $Gz" "$x0 $y0 $z0"]]set Iz [expr ($area*$distance*$distance)]lappend IZ $Izset distancex [vecmag [vecsub "$Gx $Gy $Gz" "$Gx $y0 $z0"]]set Ix [expr ($area*$distancex*$distancex)]lappend IX $Ixset distancey [vecmag [vecsub "$Gx $Gy $Gz" "$x0 $Gy $z0"]]set Iy [expr ($area*$distancey*$distancey)]lappend IY $Iy}set polarinertiamoment [sum $IZ]set inertiamomentx [sum $IX]set inertiamomenty [sum $IY]puts "截面极惯性矩为：$polarinertiamoment"puts "截面x惯性矩为：$inertiamomentx"puts "截面y惯性矩为：$inertiamomenty"
至此插件已全部完成。
04
—
单元尺寸对求解精度的影响
