基于OpenSEES虚功原理的钢结构优化

【小记】这篇是指导华南理工大学毕业设计的一篇毕业论文，今年指导三位华南理工大学的学生的毕业论文研究，三个研究都是围绕优化算法在工程上的应用展开，通过借助有限元计算软件OPENSEES及最新的OPENSEESPY （PYTHON）对结构优化算法进行实现，可用于将来的结构设计当中。每次的毕业设计论文，我都鼓励学生学习编程，通地编程的方法解决工程问题。

虚功原理的应用
虚功原理是分析力学部分的一个很重要的原理，其利用虚位移和理想约束的概念来讨论力学系统处于平衡的条件。其内容可以基本概括为：受有理想约束的力学系统处于平衡的充要条件是作用在系统上的所有主动力在任意虚位移下所作的虚功之和为零。同时，虚功原理建立了构件尺寸与位移的关系，由虚功原理得位移计算的一般公式为：
  
该公式可以理解为结构中每个单元均对施加单位力处的节点的位移做出了贡献，现定义敏感指数SI，其表达式为：
 
V为结构中每个构件的体积，△/Vmax为贡献的位移值除以体积后的最大值。同时，SI也是ETABS与SAP2000中“显示虚功图”中数值（例如下图中的结构所示）的原理。
 
由于敏感值是相对值，我们可以通过敏感值的大小对构件截面进行调整，以达到结构优化的目的。对于桁架结构，由于只存在轴力，故我们可以只看公式中的第一项。

 
当构件的SI较大时，提高该截面面积，便可以有效降低该构件带来的△值，同时，SI越大，位移降低的效果也就越明显。
该原理比较容易理解。在进行结构优化时，需要多次迭代计算，因此采用无界面的OPENSEES是不错的选择。同时在网上也并没有将OPENSEES与虚功原理相结合的结构算例，因此本文将通过一个简单的桁架实例来介绍该方法的具体应用。

具体优化算例
实例：对一总体积恒定的钢桁架结构进行结构优化，使桁架在荷载作用下（荷载布置下如图所示）的变形最小，即结构刚度最大。图中F=50000N，所有构件的截面初始尺寸均为10000 mm2，桁架结构总体积为824558cm3。单位力大小设为1000N，施加在节点11处，方向竖直向下。
 

优化流程

（1）对该实例我们采用OPENSEES进行分析。在分析之前先借助ETABS建立模型并导出sk2文件以供ETO程序生成TCL代码。ETABS中的模型如下图所示。
 

（2）导入ETO生成OPENSEES的分析命令流文件，导入ETO以后处到桁架的构件编号如图所示。在OPENSEES中我们需要记录并输出各节点的竖向位移，与各构件的内力。由于是桁架结构，输出的内力结果中只有轴力为非0项。程序求解所有节点的竖向位移的最大值，作为优化结果的评价标准，最大位移越小，则优化结果越好。 
 


（3）OpenSEES的命令流如下所示，命令流有两个部分，即实际荷载作用下的结果，单位力作用下的结果，每一部分又包括3个板块：
第一板块：节点、材料与支座条件（分析过程中不发生改变）
第二板块：截面（分析过程中不断改变），采用引入文件  source element.tcl
第三板块：分析与结果提取，需要提取节点竖向位移与各构件内力
同时，
recorder Node -file $dataDir/Disp/node1.out -time -node 1 -dof 3 disp
……
recorder Element -file $dataDir/Force/element1.out -time -ele 1 localForce
……
不同的节点或构件分析存储的文件名不同， 分别是 node1~node14.out与element.1~element25.out。
基本TCL命令流如下：
wipe
model basic -ndm 3 -ndf 3
source node.tcl
fix 1 0 1 0 0 0 0;
……
uniaxialMaterial Elastic 1 1.999E+005
source element1.tcl
## Load Case = Live
pattern Plain 1 Linear {
load 1 0.000E+000 0.000E+000 -5.000E+004 0.000E+000 0.000E+000 0.000E+000
……
}
puts "recorder1"
set dataDir Liveoutput
file mkdir $dataDir
file mkdir $dataDir/Disp
file mkdir $dataDir/Force
recorder Node -file $dataDir/Disp/node1.out -time -node 1 -dof 3 disp
……
recorder Element -file $dataDir/Force/element1.out -time -ele 1 localForce
……
constraints Plain
numberer Plain
system BandGeneral
test EnergyIncr 1.0e-6 200
algorithm Newton
integrator LoadControl 1.000E-002
analysis Static
analyze 100
……（与上一部分相同）
## Load Case = Virtual
pattern Plain 1 Linear {
load 11 0.000E+000 0.000E+000 -1.000E+003 0.000E+000 0.000E+000 0.000E+000
……
}
set dataDir Virtualoutput
file mkdir $dataDir
file mkdir $dataDir/Disp
file mkdir $dataDir/Force
……（与上一部分相同）
wipe

（4）桁架截面在优化过程中是不断变化的，所以element.tcl文件需要不断的修改，如其中一次迭代计算的截面文件element5.tcl如下所示。
element truss 1 2 1 12306.416996915492 1
……
element truss 25 6 14 11441.348183424036 1

（5）OPENSEES文件基本编写完成，现需编写Python文件。Python编码流程大致如下图所示。每次计算分析所调用的elment.tcl文件都会储存在各自的子文件夹下，以便查看优化的最终解。
 

（6）接下来讲解Python中结构体积重分配的具体原理原理。敏感值SI是相对值，其值范围在0至100之间，故以50作为截面调整的分水岭，即杆件的SI大于50时，增大截面面积；SI小于50时，减小截面面积。调整的幅度取决于SI—50的绝对值的大小。具体的公式为：
 
其中A’为优化后的面积，A为优化前的面积。且为保证结构总体积不变，求得优化前后结构的体积差，将相差的体积按SI的大小对构件进行分配。且规定截面最小面积不得小于10cm2，最大不得大于250cm2。
这里需要对截面的最小面积进行严格控制。截面面积过小可能会导致某些杆件的刚度过小，从而使某些节点的竖向位移大于施加单位力的节点（节点11），这会导致控制节点的选取与实际不符，优化效果大打折扣。

（7）运行程序，查看前50次的运算分析结果，结构最大节点位移（节点11）的位移优化曲线如图所示。
 
由曲线可知，随着不断的优化计算，结构最大节点位移不断下降，结构刚度越来越大，最大节点的位移由3.08下降到了2.09左右，下降幅度在33%左右。且最大节点位移始终在节点11处。

（8）结构优化前各截面面积均为100cm2，优化后各构件截面面积如下图所示。
 

（9）采用ETABS对比分析：
优化前结果：节点11竖向位移为3.08mm。
 
优化后结果：节点11竖向位移为2.09mm。
 

（10）优化后的位移云图：采用PYTHON进行对OPENSEES的后处理动画